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Abstract 

Dramatic progress has been made over the last decade in the numerical study of quantum 
chromodynamics (QCD) through the use of improved formulations of QCD on the lattice 
(improved actions), the development of new algorithms and the rapid increase in comput- 
ing power available to lattice gauge theorists. In this article we describe simulations of full 
QCD using the improved staggered quark formalism, "asqtad" fermions. These simula- 
tions were carried out with two degenerate flavors of light quarks (up and down) and with 
one heavier flavor, the strange quark. Several light quark masses, down to about 3 times the 
physical light quark mass, and six lattice spacings have been used. These enable controlled 
continuum and chiral extrapolations of many low energy QCD observables. We review the 
improved staggered formalism, emphasizing both advantages and drawbacks. In particu- 
lar, we review the procedure for removing unwanted staggered species in the continuum 
limit. We then describe the asqtad lattice ensembles created by the MILC Collaboration. 
All MILC lattice ensembles are publicly available, and they have been used extensively by 
a number of lattice gauge theory groups. We review physics results obtained with them, 
and discuss the impact of these results on phenomenology. Topics include the heavy quark 
potential, spectrum of light hadrons, quark masses, decay constant of light and heavy-light 
pseudoscalar mesons, semileptonic form factors, nucleon structure, scattering lengths and 
more. We conclude with a brief look at highly promising future prospects. 

PACS numbers: 1 2.38. Gc, 11. 15. Ha 
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I. INTRODUCTION 



The standard model of high energy physics encompasses our current knowledge of the funda- 
mental interactions of subatomic physics. It consists of two quantum field theories: the Weinberg- 
Salam theory of electromagnetic and weak interactions, and quantum chromodynamics (QCD), 
the theory of the strong interactions. The standard model has been enormously successful in ex- 
plaining a wealth of data produced in accelerator and cosmic ray experiments over the past thirty 
years. Our knowledge of it is incomplete, however, because it has been difficult to extract many of 
the most interesting predictions of QCD: those that depend on the strong coupling regime of the 
theory and therefore require nonperturbative calculations. 

At present, the only means of carrying out nonperturbative QCD calculations from first princi- 
ples and with controlled errors is through large-scale numerical simulations within the framework 
of lattice gauge theory. These simulations are needed to obtain a quantitative understanding of the 
physical phenomena controlled by the strong interactions, such as the masses, widths, and scatter- 
ing lengths of the light hadrons, and to make possible the determination of the weak interaction 
Cabibbo-Kobayashi-Maskawa (CKM) matrix elements from experiment. A central objective of 
the experimental program in high-energy physics, and of lattice QCD simulations, is to determine 
the range of validity of the standard model, and to search for new physics beyond it. Thus, QCD 
simulations play an important role in efforts to obtain a deeper understanding of the fundamental 
laws of physics. 

Major progress has been made in the numerical study of QCD over the last decade through 
the use of improved formulations of QCD on the lattice, the development of new algorithms, and 
the increase in computing power available to lattice gauge theorists. The lattice formulation of 
QCD is not merely a numerical approximation to the continuum formulation. The lattice regular- 
ization is every bit as valid as any of the popular continuum regularizations, and has the distinct 
advantage of being nonperturbative. The lattice spacing a establishes a momentum cutoff n/a that 
removes ultraviolet divergences. Standard renormalization methods apply, and in the perturbative 
regime they allow a straightforward conversion of lattice results to any of the standard continuum 
regularization schemes. 

There are several formulations of the lattice QCD Lagrangian in current widespread use. The 
gauge field action can be constructed with varying degrees of improvement that are designed to 
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reduce cutoff effects at nonzer o lattic e spacing. The quark action can be formulated using Wil- 
son's original method (iWilsonl 1 19741) with modern improvements (LS heikholesl ami and Wohlerti 



1985) or with the twisted mass flFrezzotti et al. 



variants (|Morningstar and Peardon 



2004 



gered fermion formulation (|Banks et al. 



2000 



Zanotti et al. 



1976 



1977 



2001 



Frezzotti and Rossil . 



20041) or other 



200 2J), with the Kogut-Susskind or stag 



Kogut and Susskin 



1975 



Suss kind. 



19771) 



with improvem ents, and with the more recently imp l emented chira 



wall fermions (Furman and Shamiii 



([Narayanan and Neubergerl . 



1995 



1995 



Ka planl 



992; 



Shamir , 



methods that include domain- 



1993b and overlap fermions 



Neubergerl Il998b|) . Other improve ments also in production 



use are Wilson quar ks with HYP smearing to reduce lattic e artifacts (IHasenfratz et all 



2007 



Schaefer et al. . 



20071) . or to approximate good chiral behavior (IGattringeii 120011) . 



In this article, we review a ten-y ear research program founded on a particular improvement of 



staggered fermions called "asqtad" ([B ernard et al 



1999 



Lepage! 



1998 



Orginos and Toussaint , 



1999; 



2000 



Blum 



Orginos et al. 



et al. 



1997 



Lagae and Sinclair!. 



1999|) (named for its O (a 2 ) level 



of improvement and its inclusion of a "tadpole" renormalization). Over this time, the MILC Col- 
laboration has created significant library of gauge field configuration ensembles with the full com- 
plement of the light sea quarks u, d, and s. The masses of the u and d quarks have been taken 
to be equal, which has a negligible effect (< 1%) on isospin-averaged quantities. In planning the 
parameters of these ensembles, an attempt has been made to address the three primary sources of 
systematic errors in lattice QCD calculations: the chiral and continuum extrapolations and finite 
size effects. It is straightforward to perform simulations with the mass of the s quark close to its 
physical value, and in most of the ensembles that has been done. However, up to now it has been 
too computationally expensive to perform simulations at the physical mass of the u and d quarks. 
Instead, ensembles have been generated with a range of light quark masses in order to perform 
extrapolations to the chiral (physical value of the u and d quark mass) limit guided by chiral per- 
turbation theory. Simulations have been performed with six values of the lattice spacing in order 
to enable controlled extrapolations to the continuum (zero lattice spacing) limit, and in almost all 
cases the physical size of the box in which the simulations have been carried out has been taken 
to be more than four times the Compton wavelength of the pion in order to minimize finite size 
effects. Finally, because SU(3) chiral perturbation theory converges rather slowly for the s quark 
mass close to its physical value, a number of ensembles have been generated with lighter than 
physical s quark masses to improve the chiral extrapolation. These ensembles are publicly avail- 
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able, and have been used by a number of research groups, including our own, to calculate a wide 
variety of hadronic quantities ranging from chiral properties of light mesons to hadronic parton 
distributions to semileptonic decays of mesons with a charm or bottom quark to the spectroscopy 
of heavy quarkonium. 

The asqtad improved staggered fermion approach has enjoyed considerable success. Its com- 
paratively high degree of improvement and its relatively low computational cost enabled a broad 
set of full QCD phenomenological calculations earlier than was possible with other fermion meth- 
ods. In Fi g. Q] we illustrate th e dramatic effects of including sea quarks in a variety of physical 
quantities (|Davies et all 120041) . Computations with asqtad sea quarks are able to account for a 
wide variety of known decay constan ts, some hadronic m asses, and several quarkonium mass split- 



tings to a precision of a few p ercent (|Davies et al. . 



2004). Their predi c tions fo r a few heavy-light 



leptonic (|Aubin et all I2005al ) and semileptonic decays (jAubin et all l2005bT) have been experi 



menta lly confirmed. They pr ovide values for the strong fine structure constant (%<■ (IDavies et al. 



2008), the cha rm quark mass (IDavies et al. 



20091). the CKM matrix elements | V us | ([Bernard et al 



2007ej), \V cb \ (Bernard et all l2009ah and \V ub \ (IBailev et all \200% . and the D + and D s leptonic 



decay constants (|Follana et all |2008|) that are competitive with the most accurate determinations 
to date. 

In Sec. HH we begin with a brief review of lattice gauge theory, discussing gauge field and 
fermion field formulations and numerical simulation methods. We end Sec. HI] with an overview of 
the asqtad and the more recent HISQ fermion formulations. 

Section Un] first discusses the inclusion of staggered discretization errors in chiral perturbation 
theory, resulting in "staggered chiral perturbation theory" (SXPT). The application to the light 
pseudoscalar meson sector is described in detail; the applications to heavy-light mesons and to a 
mixed-action theory (with chiral valence quarks and staggered sea quarks) are treated more briefly. 
We then turn attention to the procedure we use to deal with the extra species that occur for stag- 
gered fermions. Each staggered field (each flavor of quark) normally gives rise to four species 
in the continuum limit. The additional degree of freedom is called "taste." To obtain the correct 
counting of sea quarks it is necessary to take the fourth-root of the fermion determinant. This 
rooting procedure has been shown to produce a theory that is nonlocal on the lattice, leading to 
the legitimate question of whether the nonlocality persists as the lattice spacing goes to zero. Such 
nonlocality would spoil the continuum limit, giving a theory inequivalent to QCD. In recent years, 
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FIG. 1 Comparison of the ratio of lattice QCD and experimental values for several observables, where the 
lattice QCD calculations are done in the quenched approx imation (left) and w ith 2 + 1 flavors of asqtad sea 



quarks (right). This is an updated version of a figure from 



Davies et al. 



(2004). 



however, there has been a considerable amount of work on this issue, and there is now a substantial 
body of theoretical and computational evidence that the fourth-root methodology is indeed correct. 
We discuss some of that work in detail in Sec. [nil and also explain how to take rooting into account 
properly in the chiral effective theory. 

In SecUVj we list the ensembles of publicly available asqtad gauge configurations, and describe 
tests of their intended properties, including the determination of the lattice scale and the topological 
susceptibility. In the following sections, we review physics results obtained with them. In Sec.lVl 
we review the spectroscopy of light hadrons other than the pseudoscalar mesons, including vector 
and scalar mesons and baryons. Section[VI]is devoted to properties of the pseudoscalar mesons, in- 
cluding masses, decay constants and Gasser-Leutwyler low energy constants. We turn in Secs. lVIIl 
and IVIIII to the masses and decays of mesons containing one heavy (charm or bottom) quark and 
one light antiquark. Section rvTll treats masses and leptonic decays; Sec. I Villi semileptonic decays. 

In Sec. [DO we review a variety of other calculations, including the determination of the strong 
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coupling 0C5, quarkonium spectroscopy, the spectroscopy of baryons containing one or two heavy 
quarks, — and Bq — Bq mixing, the muon anomalous magnetic moment, and quark and gluon 
propagators. 

Finally, in Sec.|Xl we discuss further improvements under way or under consideration, including 
the incorporation of electromagnetic effects and the implementation of the HISQ action, and briefly 
comment on future prospects for the field. 

We do not review applica tions of the asqtad form ulation to QCD thermodynamics. A recent 



article by DeTar and Heller (IDeTar and Hellerl . 



20091) contains a review of high temperature and 



nonzero density results, including those obtained using the asqtad fermion action. 

II. FERMIONS ON THE LATTICE: IMPROVED STAGGERED FORMALISM 
A. Brief introduction to lattice gauge theory 

1. Basic setup 

Euclidean, i.e., imaginary time, field theories can be regulated by formulating them on a space- 
time lattice, with the lattice points, called sites, separated by the lattice spacing a. This introduces 
an ultraviolet cutoff n/a on any momentum component. Matter fields then reside only on the 
lattice sites, while the gauge fields are associated with the links joining neighboring sites. The 
gauge fields are represented by gauge group elements UJx) on the links, which represent parallel 
transporters from site x to the neighboring site x + ajj, where ft is the unit vector in the direction /j, 
with /j = 1 , . . . , d for a J-dimensional lattice: 

Un(x) = iPexpjigjf dy v A v (y)^j =exp!^iga A^x + afr/l) + ^dfa^x + afr/Z) + . . . J 
= l + iagA tl (x+ajj/2) + ... . (1) 

Under gauge transformations V(x), restricted to the sites of the lattice, the gauge links transform 

as 

U p (x) -> V(x)U,(x)V\x + aju) . (2) 

The traces of products of gauge links around closed loops on the lattice, so-called Wilson loops, are 
then gauge invariant. The gauge action can be built from the sum over the lattice of combinations 
of small Wilson loops with coefficients adjusted such that in the continuum limit, a — > 0, it reduces 
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to f d d x\TrF^ y up to terms of (a 2 ). The simplest gauge action, the original action introduced by 



Wilson ( 



19741) . consists of a sum over plaquettes (lxl Wilson loops) 



5 G = ^£ReTr(l-t/ pZ ), (3) 

pi 

where |3 = 2N/g 2 , for gauge group SU(iV), with g 2 the bare coupling constant. 

Fermions, in Euclidean space, are represented by Grassmann fields \\r x and which in the 
lattice formulation reside on the sites of the lattice. A generic fermion action can be written as 

Sf = Y^xMp-^yVy , (4) 

x,y 

where the fermion matrix Mf- XJ is some lattice discretization of the continuum Dirac operator 
D + m. Details of lattice fermion actions are described below. 
The lattice gauge theory partition function is then given by 

Z(P) = / JJdU^x) nt^MV*] exp{-5 G - a 4 S F } , (5) 

X,\l X 

where dU^(x) is the invariant SU(A0 Haar measure and d\\f x d\\f x indicate integration over the Grass- 
mann fields. 

Since Sf is quadratic in the fermion fields, the integration over the Grassmann fields can be 
carried out, leading to (up to a trivial overall factor) 

Z(P) = f n^WdetM F exp{-5 G } = ! n^(*)exp{-S e// } , (6) 

with S e ff = So — TrlogM/7. 

The expectation value of some observable O is given by 

(O) = ^J]]dU^x)]][dy x dy x }Oexp{-S G -a 4 S F } 

= ^ I Y[dU,{x)OdtiM F ^ V {-S G } = / Y\dU,(x)Oex V {-S eff } . (7) 

If the observable O involves fermion fields \\f x and \jjfy then, in the second line of Eq. © each 
pair is replaced by Mp l xy in all possible combinations with the appropriate minus signs for Wick 
contractions of fermion fields. 
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2. Improved action 



As mentioned before Eq. ©, the typical gauge action on the lattices reduces to the continuum 
action up to terms of o(a 2 ). These terms lead to 0(a 2 ) deviations from the continuum result 
of physical observables computed at finite lattice spacing. These (a 2 ) effects can be reduced 
by using an improved gauge acti on (together with impro ved operators, where necessary) in an 



1983) 



improvement program initiated by iSymanzikl (| 19801 . 

For the gauge action, the improvement can be achieved by adding 2x1 (planar) rectangle 
(labeled "rt") and generalized 3-d all 1 x 1 x 1 parallelogram (labeled "pg") Wilson loop terms (see 
Fig. [2]) to the Wilson action, Eq. Q, w ith coefficients computed, at one-loop order in perturbation 



theory, by 



Luscher and Weisz (1985a 



bl 



S w = | \ J> p/ ReTr(l - U pl ) +£c rt ReTr(l - U rt ) +£c pg ReTr(l - U pg ) \ 
7V [pi rt P g J 



(8) 



he coe fficients, c, = c- "* +4%0Lqc^ at one loop, can be found in Table 1 of 



Luscher and Weisz 



(1985a). 



A 



a) ^ b) V- c) V 

FIG. 2 Luscher- Weisz action Wilson loops: a) standard plaquette, b) 2 x 1 rectangle and c) 1 x 1 x 1 

parallelogram 



Bare lattice perturbation theory resul ts generally converge slowly bu t can be improved by using 
tadpole-improved perturbation theory (|Lepage and MackenzieL 1 19931) . This starts with using a 
more continuum-like gauge link U M — > = Uq U^. The so-called tadpole factor uq is determined 
in numerical simulations either as the expectation value of U p in Landau gauge or, more commonly, 
from the expectation value of the average plaquette 



Uq 



(iReTrtV) 1 / 4 . 



(9) 



The Luscher- Weisz action can now be tadpole improved by explicitly pulling a u Q 1 factor out of 
each link and replacing oco in the one-loop perturbative coefficients c ; with a nonperturbatively 
renormalized coupling a s defined, for gauge group SU(3), in terms of the measured lattice value 
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of uq by 



a, = -1.3036 15 log wo 



(10) 



where the proportionality factor is determined by the one-loop expression for log wo- Defining 
3/,w = MQ 4 3cn/, sinc e U p i involves the product of four links, the improved action can be written as 
(lAlford et all 



1995) 



Votm r; "\ v [l+0.4805aj D v 0033250, 
2^ReTr(l - t/ p/ ) - ^ ^ ReTr ( 1 ~ ~L ^ ^^C 1 ~ U ps> 



Pi 



pa 



(ID 



Since higher perturbative orders in the coefficients are neglected, the one-loop improved Liischer- 
Weisz action, Eq. CCD, leads to remaining lattice artifacts of O (a 2 a 2 ). Sometimes, only a tree-level 
improved action without the terms proportional to a, in Eq. (fTT|) is used, leading to lattice artifacts 
of O (oc ? a 2 ) . Since the parallelogram terms are then absent such simulations are somewhat faster. It 
should be noted that Eq. (ITTT) does not include the one-loop contributions from dynamical fermions, 
which were unknown at the time the MILC collaboration started the 2 + 1 -flavor simulations re- 
viewed in this article. Therefore, for those simulations, the leading lattice artifacts in the gauge 
sector are O (oc 5 a 2 ) as in th e fermion sector, described later. The one-loop fermion contribution has 



recently been computed by 



Hao. von Hippel. Horgan. Mason, and Trottierl (|2007D . 



B. Fermions on the lattice 

1. The doubling problem 

Putting fermions on a lattice, one replaces the covariant derivative in the continuum fermion 
action with a covariant (central) difference 

Snaive = £?W |EYa/ V a/VW +m\\f(x) J , (12) 

where 

VM X ) = Ya (^MV^ + aA) -tf* • (13) 
The inverse propagator in momentum space derived from the action Eq. (fT2l) in the free case, with 
all link fields U p = 1, is 

aS~ l (ap) = z'^y jU sin(ap jU ) +am . (14) 
13 



In the massless case, this inverse propagator not only vanishes when p = 0, but also when p^ = 
or p ll = Tl/a for each /j = 1, ... ,4, i.e., on all 16 corners of the Brillouin zone in d = 4 dimensions. 
Thus, when we try to put one fermion on the lattice we actually get 16 in the continuum limit. This 
is the infamous doubling problem of lattice fermions. 



2. Wilson fermions 

This doubling problem was recognized by Wilson when he first formulated lattice gauge theo- 
ries. He also prop osed a solution : adding an irrelevant term — a term that vanishes in the contin- 



uum limit, a — > (|Wilson , 



1975) 



(XV 

S W = Snaive ~ "y Y, A M X ) = , (15) 

Z X fl 

where r is a free parameter, usually set to r = 1, and the Laplacian is 

A (U \)/(jc) = — j {u^ (x) \\f(x + afj) + Uj, (x - api) \|/(x - ap) - 2\|/(x)^ . (16) 
The free inverse propagator now is 

aS~ l (ap) = i^y^smfap/j) +am — (cos(ap (U ) — l) . (17) 

The doublers, with n momentum components p^ = %/a, now attain masses m + 2nr/a, and only 
one fermion, with p ps 0, remains light. 

We note that the Wilson Dirac operator is ys-Hermitian, 

Dl(m)=j 5 D w (m)j 5 . (18) 

Thus detD^m) = det Dw(m), implying that two flavors — and by extension any even number of 
flavors of Wilson fermions — lead to a manifestly positive (semi-) definite fermion determinant, 
det[Dl(m)D w (m)}. 

The price for eliminating the doubling problem in this Wilson fermion approach is that the 
action Eq. (fT5l) violates the chiral symmetry 8\|/ = iujsy, 8\j/ = j'caj/ys of massless fermions (with 
a an infinitesimal parameter). As a consequence, the massless limit of fermions is no longer 
protected — the mass gets an additive renormalization; to get massless quarks requires a fine 
tuning of the bare mass parameter. 

14 



According to the usual, renormalization group based universality arguments, the chiral symme- 
try, broken at finite lattice spacing only by an irrelevant, dimension-five operator, will be recovered 
in the continuum limit after fine tuning of the bare mass parameter. But the explicit violation of 
chiral symmetry allows the generation of other contributions to dimension-five operators which are 
suppressed by only one power of the lattice spacing a. The lattice artifacts for Wilson fermions are 
therefore of (a), rather than O (a 2 ) as in the pure gauge sector. 

Besides vj/(jc)A\|/(x), with A = L^A^, there is a second dimension-five (chiral symmetry break- 
ing) operator 

Ssw = ^cswY,V( x ) a t"?i*>( x )V( x ) , (19) 



where f^ipc) is a lattice representation of the field strength tensor F^ix), and c^ v = jpY^Yv]- 
In clusion of Eq. (tT9b into the fermion a ction, with properly adjusted coefficient csw, was proposed 
by lSheikholeslami and Wohlertl (11985b to eliminate the O (a) effects of the Wilson fermion action. 
Since J„ v (x) on the lattice is usually represented by a "clover leaf" pattern of open plaquettes, the 
action including the term Eq. (fT9l is commonly referred to as the clover action. 

The appropriate coefficient csw of the clover term, Eq. (fT9l , can be computed in perturbation 



1996. 



theory (ILuscher and Weisz , 



1996 



Wohlert, 



19871) . or even better, nonperturbatively (|Luscher et all 
19971) - truly reducing the remaining lattice effects from (a) to (a 2 ) . 



Another problem with Wilson fermions is that, because of the additive mass renormalization, 
the fermion determinant det£%(m) is not positive definite even for putative positive quark mass. 
Configurations with detDiy(m) ~ can occur, called exceptional configurations, which can slow 



down numerical simulations considera bly. A formulation that removes such exceptiona 



urations, introduced by Frezzotti et al. ([Frezzotti et al. 



2000, 



2001; 



Frezzotti and Rossi. 



confi g 



2004|) is 



called "twisted-mass QCD". For two flavors one considers the Dirac operator 



D t 



D + m + Z//Y5T3 



(20) 



where the isospin generator X3 acts in flavor space. In the continuum, the twisted-mass Dirac op- 



erator is equivalent to a usual Dirac operator with mass >Jm 2 + fj 2 . On the lattice, however, with 
D replaced by the (massless) Wilson Dirac operator £V(0) of Eq. (fT5l) . the twisted-mass term 
ensures a positive-definite two-flavor determinant, det[D^(m)Z)^(m) +/J 2 ] > 0. An added ben- 
efit of the twisted-mass (Wilson) fermion formulation is, that at max imal twist tana = u/m, th e 



twisted-mass Wilson Dirac operator is automatically O (a 2 ) improved (Frezzotti and Rossi . 



2004). 
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Unfortunately, the real part of the mass m still receives an additive renormalization so that achiev- 
ing maximal twist requires a fine tuning. Furthermore, at finite lattice spacing, isospin symmetry 
is broken, making the 7t° mass different from the mass of the n ± . 



3. Staggered fermions 



Another way of dealing with the doubling problem, a 



the st aggered fermion formalism (|Banks et all 
1977b . One introduces a new fermion field by 



197 



1977 : 



l eviating though not eliminating it, is 



Kogut and Susskindl 



1975 



Susskind, 



VM = T x x(x) , Vjjr(x) = x(x)T 



t 

X > 



with 



Using rjr x = 1 and 



the naive fermion action, Eq. (fT2|) . can be written as 



r — J. x i/ a ) J. x 2/a) Jib/a) Axa/o) 
1 x — J i 12 S3 14 



Sks = \ v ^ XW + ™XM \ = X (DKS + m) 1 , 



(21) 



(22) 



(23) 



(24) 



where matrix multiplication is implied in the final expression. Here, the four Dirac components 
decouple from each other, and the fermion field %(x) can be restricted to a single component, 
thereby reducing the doubling by a factor of four, from sixteen to four. It is, in principle, possible 
to interpret these four remaining degrees of freedom as physical flavor (u, d, s, c), but, in order to 
give different masses to the flavors, one mus t introduce general mass terms coupling nearby sites 



(IGockelei . 



1984 



Golterman and S mit. 



198 4). That approach then leads to a variety of practical 



problems including complex determinants and the necessity of fine tuning. 

Instead, we follow modern usage and refer to the quantum number labeling the four remaining 
fermion species as "taste," which, unlike flavor, is an unwanted degree of freedom that must be 
removed. We describe how this removal is accomplished by the so-called "fourth-root procedure" 
at the end of this section, and discuss it in more detail in Sec. lIIICl If more than one physical flavor 
is required, as is, of course, the case for simulations of QCD, one then needs to introduce a separate 
staggered field for each flavor. For example, for QCD with three light flavors, one employs three 
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staggered fields, % u , %j, and % s . 1 However, for simplicity, we consider only a single staggered field 
(one flavor) in the remainder of this section. 

The one-component fermions with action Eq. (1241) are referred to as (standard) staggered or 
Kogut-Susskind fermions. The "standard" distinguishes them from improved versions, described 
later on. 

An important discrete symmetry of the staggered ferm ion action, Eq. (1241) . is shift symmetry 



( van den Doel and Smit. 



1983 



Golterman and Smut 



1984) 



U v (x) — > U v (x + ajj) , 



(25) 



with the phase p jU (x) defined by p^x) = ( — 1)( x /j+i h i- x 4)/« Additional discrete symmetries of the 
staggered action are 90° rotations, axis inversions, and charge conjugation. In the continuum limit, 
these symmetries are expected to enlarge to a direct product of the Eu clidean Poincare group and a 



Golterman and Smiti 



1984). 



vector SU(4)v among the tastes (plus parity and charge conjugation) ( 

For massless quarks, m = 0, the staggered fermion action also has a continuous even/o dd 
U(l) e xU(l) chiral symmetry ([Kawamoto and Smit . 



1981 



Kluberg-Stern et al. 



19811 



1983b ). a 



remnant of the usual chiral symmetry for massless fermions in the continuum. The U(l) e xU(l) 
symmetry is 



%(x) -> exp{z'a e }%(jc) , %(x) -> %(x) exp{-/a } for x = even 
X(x) -> exp{ia a }x(x) , % ( x ) ^ X W exp{-z'a e } for x = odd , 



(26) 



where a e and a are the symmetry parameters, and a site x is called even or odd if ^{x^/a) is 
even or odd. The "axial pa rt" of this symmetry, a e = —a a = cc e , is known as U(l) e symmetry 



(Kawamoto and Smit . 



19811) and takes the form 



X(x) -> exp{ia e e(x)}x{x) , %{x) -> %{x) exp{/a e e(x)} with e(x) = (-l)^( x f/ a ) . (27) 

The chiral symmetry, Eq. d26l) or Eq. (1271) . protects the mass term in Eq. d24l) from additive 
renormalization, while the discrete symmetries (especially shift symmetry) are also needed to 



1 In practice, since one usually takes m u = ^ m s , the u and d fields can be simulated together, and one can use 
only two staggered fields. For clarity, we ignore this technical detail in our exposition. 
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preve nt other mass terms (coupling % and % at nearby sites) from arising (IGolterman and Smiti 
984 ). In particular, an alternative version of staggered quarks called the "Dirac-Kahler action" 



(|Becher and Jooslll982) does not have s. 



loop even when m = (IMitra and Weisz , 



lift sym metry and therefore generates a mass term at one 



1983) 



The even/odd symmetry is spontaneously broken to the diagonal vector U(l V (quark number) 
symmetry, a e = a a , with an ensuing Goldstone boson. In addition, the mass term breaks the 
U(l) f xU(l) symmetry explicitly, giving mass to the G oldstone boson, <~< m. 



The staggered Dirac operator D^s m Eq. (1241) obeys (|Smit and Vinkl 



1987) 



D 



KS 



-D KS = e%8 



(28) 



where £ is a diagonal matrix in position space with e(x) along the diagonal, and the second equality 
follows from the fact that Dks connects only even and odd sites. The fact that Dks is antihermitian 
implies that its eigenvalues are purely imaginary; the £ relation in Eq. (|28l) then tells us that the 
nonzero eigenvalues come in complex-conjugate pairs. For m > 0, which is the case of interest 
here, this ensures that the staggered determinant det(DKS + m) is strictly positive. 2 Note that 
the continuum Euclidean Dirac operator D cont is also antihermitian and obeys a corresponding 
equation 



D 



cont 



-D r 



= Ys D cont Ys , (29) 

which similarly (but now only formally) results in a positive determinant for positive quark mass. 

The one-component staggered fermion fields %(x) can be assembled into Dirac fields q(y), 
living on 2 4 hypercubes of the original lattice, labeled by y, with corne rs x = 2y + aA, where 



1982; 



GliozziL 



1982 



Kluberg-Stern et all 



1983a ). One has 



ifj = 0, 1 (|Duncan et all 

q(y)ai = ^ETOca U A (y) %(2y + aA) , q(y) ia = l£^(2y + aA) ufo) (T A )] a , (30) 



8 



A "A 

where a, i label the Dirac and taste indices, respectively, and [/^(j) is a product of the gauge links 
over some fixed path from 2y to 2y + aA. Bilinear quark opera tors, with spin structure y s = T s and 



taste structure % t = T f * are defined by (ISharpe and Patel . 



1994) 



O st = q(y)(y s ®$t)q(y) = ^Y,%(2y + aA) Ul(y) U B {y) %(2y + aB) ^tr (rl^rj) . (31) 



2 We do not expect any exact zero modes on generic configurations, even those with net topological charge. Such 
configurations will in general have only some near-zero {0(a) or smaller) eigenvalues. So, in fact, the determinant 
should be positive even for m < 0. This is different from the case of chiral fermions discussed in Sec. lII.B.4l 
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In the free case (all U fl {x) — 1) , the qu ark action in Eq. (1241) can be expressed in terms of the 



fields q{y) as (|Kluberg-Stern et al 



1983ah 



Sks = 16£<?(y) jm(/<g)/) + £ [(y„ ®l)y^ + a (y 5 ® ^ 5 ) Aj j $ (y) , (32) 

where / is the identity matrix, the factor of 16 arises from the fact that there are 1/16 as many y 
points as x points, and V„ and A^ are the free-field versions of Eqs. (fT3l) and (fT6l) . but acting on the 
doubled (y) lattice: 

= ^if(y+2w)-f(y-2w)} , 
A »f(y) = ^2 \f(y + 2a ® - 2 f(y) + f(y - 2a M ■ w 

These derivatives go to d(jf(y) and dj l f(y), respectively, in the continuum limit. In the interacting 
case there is another dimension-five, (a), term, involving the field- strength tensor in addition 

here are also higher contributions of o(a 2 ) starting at dimension six 



to the A t , term in Eq. (132 



(|Kluberg- Stern et all 



1983a). 



In the V„ (first derivative) kinetic energy term of Eq. (1321) . the even/odd U(l) e xU(l) symmetry 
is enlarged to a full continuous chiral symmetry, U(4)^xU(4)^, acting on the taste indices of the 
right and left fields, qit(y) = \ {\ + J5)q(y) and <?z,(y) = \ {\ —J5)q(y)- The mass term breaks this 
down to an SU(4)y vector taste symmetry (plus the U(l)v of quark number). On the other hand, 
because of the explicit taste matrices, the second derivative term in Eq. (l32l) breaks the full chiral 
symmetry to the U(l) e xU(l) symmetry (plus the discrete staggered symmetries). Because these 
are all symmetries of the original staggered action, they remain symmetries in the taste basis, even 
when the additional terms that appear in Eq. (1321) in the interacting case are taken into account. 

The key point is that, in the interacting theory, one can split the staggered Dirac operator in the 
taste basis as: 

D KS = D®I+aA, (34) 

where / is here the (4 x 4) identity matrix in taste space, and A is the taste-violating (traceless) part, 
with minimum dimension five. One expects the SU(4)y vector taste symmetry to be restored in the 
continuum limit because A should be irrelevant in the renormalization-group sense. 



In the free case, the shift symmetry, Eq. (|25T) . takes the form for the Dirac fields q(y) (|Luo . 
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19971) : 

l{y) -> 5((/®^ + Y5Y A i®^5)9(y) + (/®^.-Y5Y A i®^5)9(y + 2^)) , (35) 

qiy) -> 2(^)( 7 ®^-W®^)+9(y+2^)(^®^i+Y5Y«®^5))) • ( 36 ) 

As the continuum limit is approached, shifts become simply multiplication by the taste matrix 
plus higher-dimension terms involving derivatives. Thus shifts are basically discrete vector taste 
transformations, coupled with translations. 
In the taste basis, the axial U(l) e symmetry is 

qiy) -> exp {/a e (y 5 <g> 4s) } (y) , qiy) -> <?0) e*P {?'a e (y 5 <g> £5) } • (37) 

Because of the ^5, this is clearly a taste nonsinglet axial symmetry, and hence is nonanomalous. 
The anomalous axial symmetry U(1)a must be a taste-singlet: 

q(y) ^exp{ia A (y 5 ®I)}q(y) , q{y) -»• q(y)vq?{ia A (y 5 ®/)} . (38) 

Indeed, this symmetry is not an invariance of the staggered lattice action in the massless limit, 
and the symmetr y violations generate , throu gh the triangle graph, the correct axial anomaly in the 



continuum limit (Sharatchan dra et all 



119811) . 

The bilinear quark operators in Eq. (I3TT) can create (or annihilate) mesons. Therefore, for stag- 
gered quarks, each meson kind with given spin (Dirac) structure F s (e.g., T s = 75 for the pion, 
r$ = Jk for the rho, etc.) comes in sixteen varieties, labeled by the taste index t. In the contin- 
uum limit all nonsinglet mesons of a given spin are degenerate 3 - SU(4V taste symmetry connects 
them. But at nonzero lattice spacing, there is only the staggered symmetry group, the group of 
the discrete symmetries of the staggered action (shifts, 90° rotations, axis inversions, charge con- 
jugation) plus the U(l)v of quark number, which are remnants of the continuum Poincare, taste 
SU(4)y, quark number, and discrete symmetries. Meson states may be classified under the sub- 
group of the staggered symmetry group , the "staggered res t frame symmetry group," which is the 



symmetry group of the transfer matrix (|Golterman , 



1986aUbl) . The sixteen tastes of a meson with 



3 Mesons that are singlets under taste and any additional flavor symmetries need not be degenerate with the nonsinglet 
mesons, since they can have physically distinct disconnected contributions to their propagators. The most important 
example is the T|', which will get a contribution from the anomaly and have a mass in the continuum limit different 
from that of all other pseudoscalars. 
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given spin structure are not degenerate at finite lattice spacing, but are split according to irreducible 
representations of the rest frame group. In particular, only the pion with pseudoscalar taste struc- 
ture £ f = 7^ is a Goldstone boson, denoted by Tip (P stands for pseudoscalar taste), whose mass 
vanishes for massless quarks, m = 0. To leading order in the chiral expansion (see Sec. IIII.AI) the 
other tastes have masses 

ml t = ml p + a\ = 2Bm + a\ , (39) 

with B a low energy constant and d t a taste-dependent splitting that is independent of a (up to 
logarithms) for small a. The non-Goldstone pions become degenerate with the Goldstone pion 
only in the continuum l imit. The taste viola tions in the pion system are found to be larger than 



those for other hadrons dlshizuka et al. 



1994). 



Since staggered fermions have only one (spin) component per lattice site, and since they have a 
remnant chiral symmetry that insures positivity of the fermion determinant at positive quark mass, 
they are one of the cheapest fermion formulations to simulate numerically. The main drawback is 
the need to eliminate the unwanted extra tastes, using the so-called "fourth-root procedure." Each 
continuum fermion species gives a factor of detMp in the partition function, Eq. ©. Therefore, 
to reduce the contribution from four tastes to a single one, we take the fourth root of the determi- 
nant, (det M ks) ] ^ , where Mrs = Dxs + m®I, with Dks given in Eq. (l34l) . The procedure was 
first introduced in the two dimensional version o f staggered fermions (where it is a "square-root 



procedure" because there are only two tastes) by 



Marinari. Parisi. and Rebbil (|1981br) . The point 



is that the Dirac operator Dj^s (and hence Mks) should become block diagonal in taste space in 
the continuum limit because A is an irrelevant operator. The fourth-root procedure then becomes 
equivalent to replacing the D^s by its restriction to a single taste. Conversely, the nontriviality of 
the prescription arises because taste symmetry is broken at nonzero lattice spacing. In Sec. IIII.CL 
we discuss the status of this procedure and the evidence that it accomplishes the goal of producing, 
in the continuum limit, a single quark species with a local action. 



4. Chirally invariant fermions 

None of the ways of dealing with the fermion doubling problem outlined so far are entirely 
satisfactory. Wilson-type fermions explicitly break chiral symmetry, and staggered fermions have 
a remaining doubling problem, requiring the fourth-root procedure, that continues to be somewhat 



21 



controversial because of the broken taste symmetry at finite lattice spacing. 

Indeed, the chiral anomaly impl ies that no lattice action can hav e an exact flavor-singlet chira l 



symmetry (IKarsten and Smitlll98ll) . There is even a no-go theorem ([Nielsen and Ninomiya , 



1981) 



that states that the doubling can not be avoided with an ultralocal 4 and unitary fermion action. 
However, actions with a modified form of chiral symmetry on the lattice can avoid doubling while 
retaining most of the desirable features of chiral symmetry. Such actions couple arbitrarily distant 
points on the lattice but with exponentially suppressed couplings, exp{— r/r^}, where should 
be of the order the lattice spacing to ensure a local action in the continuum limit. There are three 
known ways of achieving this. 



The first go es u nder the name of "domain- w all fermions" and was developed by lKaplanl(| 19921) . 



Shamir ( 



1993), and 



Furman and Shamir! (| 19951 ). The construction of Furman and Shamir is usually 
used nowadays. One introduces an additional, fifth dimension of length L s and considers 5-d 
Wilson fermions with no gauge links in the fifth direction, and the 4-d gauge links independent of 
the fifth coordinate s, 

Sdw = £ |E ~ ^ A <") VC*.*) ~P-V(x,s+ 1) -P+y(x,s- 1) 

(40) 

where P± = j(l ± Ys) are chiral projectors and we have set r = a = 1. M, introduced here with a 
sign opposite that of the mass term for Wilson fermions (fT5l) . is often referred to as the domain- 
wall height and needs to be chosen such that < M < 2. For free fermions, M = 1 is the optimal 
choice, while in the interacting case M should be somewhat larger. The fermion fields satisfy the 
boundary condition in the fifth direction, 



P-y(x,L s ) = -m f P-y(x,0) , -1) = -m f P+y(x,L s - 1) 



(41) 



where rtif is a bare quark mass. 

For ntf = 0, the domain- wall action, Eq. (|40l) . has 4-d chiral modes bound exponentially to the 



4 We denote by "ultralocal" an action that couples only sites a finite number of lattice spacings apart. A "local" action 
is either ultralocal, or the coupling falls off exponentially with distance with a range of the order of a few lattice 
spacings, so that the action becomes local in the continuum limit. Such "local" actions are believed not to change 
the universality class in the renormalization group sense. Any other action is called "nonlocal." 
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boundaries at s — and s = L s — 1 , which are identified with the chiral modes of 4-d fermions as 



q R (x) = P+y(x,L s - 1) , q L {x) = P-V(x,0) , q R (x) = \j/(x,L, - 1)P_ , q L (x) = y{x,0)P+ • 

(42) 

When L s — > °o the chiral modes become exact zero modes, the left and right handed modes q L 
and q R do not interact for nif = 0, and the domain- wall action has a chiral symmetry. At finite L s 
the chiral symmetry is slightly broken. Often L s = (10 — 20) is large enough to keep the chiral 
symmetry breaking negligibly small. The computational cost of domain-wall fermions is roughly 
a factor of L s larger than that for Wilson-type fermions. 



Naravanan and Neubergei 


(1995 


); 


Neuberger 


( 1998b 


). The overlap Dirac operator for massless 


fermions can be written as 


(Neubergei 




1998b 





aD 0V = M[l+y 5 ®(j 5 Dw(-M))} 



(43) 



where Dw(— M) is the usual Wilson Dirac operator with negative mass m = —M, and again < 
M < 2 should be used. &(X) is the matrix sign function, for a Hermitian matrix X, that can be 
defined as 

&{X) = -J= . (44) 



X 2 

Using the fact that ® 2 (X) = 1, i t is easy to see that the Neub erger Dirac operator satisfies the 



so-called Ginsparg-Wilson relation (|Ginsparg and Wilsonl 



1982), 



{j5,D m ,} = aD ov y 5 RD ov , 
with R = l/M, or equivalently, when the inverse of D ov is well defined, 

{y5,D~}} =ay 5 R. 



(45) 



(46) 



In the continuum, chiral symmetry implies that the massless fermion propagator anticommutes 
with Y5. The massless overlap propagator violates this only by a local term that vanishes in the 
continuum limit. Ginsparg and Wilson argued that this is the mildest violation of the continuum 
chiral symmetry on the lattice possible. In fact, any Dirac operator sa tisfying the Gin sparg-Wilson 
relation (1431 ) has a modified chiral symmetry at finite lattice spacing (|Luscherlll998|) . 



8y = iay 5 (l - y , 6y = icrifr (l - ^d) y 5 . 



(47) 
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or 



8\|/ = jay 5 ( 1 - J y = iay 5 \\f , 8\j/ = za\j/y 5 , (48) 

with Y5 = Y5 (l — jjD) satisfying = 75 and, using the G-W relation, Eq. (l45l) . y| = 1 . 

The close connection between domain-wall and overlap fermions can be made more ex- 
plicit by integrating out the "bulk fermions", whi ch hav e mass e s of the order of th e cutof f 
1/q, from the domain- wall a c tion. Eg. d40l). see 



Kikukawa and Noguchi 



41999b; 



Neubergei 



Boricil d2000h 



Edwards and Helled d200lh : 



1998a ). In the limit L s —> °°, one ends up with the 



overlap Dirac operator, but with the Hermitian Wilson kernel Hw = JsDw in Eq. (l43l) replaced by 
a more complicated Hermitian kernel, 

1 1 



H 7 



-H 



w 



H 



w 



(49) 



\+2a 5 H w j 5 ,v \+2a 5 H w y 5 

Here we denote the lattice spacing in the fifth direction by a$. It is usually chosen to be the same 
as the 4-d lattice spacing, as = a, which, in turn, is usually set to 1. From Eq. (|49l) we see that 
domain-wall fermions in the limit L s — > °°, followed by the limit a$ — > become identical to overlap 
fermions with the standard Neuberger Dirac operator. 

The difficulty with numerical simulations using overlap fermions is the evaluation of the sign 
function ©(Hw) of the Hermitian W ilson Dirac op erator Hw = Ys^w in Eq. (1431) . This can be 
done with a Lanczos-type algorithm (|Boricil . Il999|) . Alternatively, &(Hw) can be represented as 
a polynomial, or, mor e efficiently, as a rat ional function that can be rewritten as a sum over poles 



([Edwards et all 



1999 



Neuberger, 



Zolotarev, first given in 



1998 a|), with th e optimal approximation, using a theorem of 



van den Eshof et al. 



(l2002h . 



®(Hw)=H w LjClj W 2j =H w 



Ck 



All d^s are positive, and the necess ary inversions with the 



multishift conjugate gradient inverter (|Frommer et al. 



1995 : 



\H w +d k 



sparse matrix H 



Jegerlehner 



1996 



(50) 



are d one using a 



1998) 



Finally, two versions of fermions that satisfy the Gi nsparg- Wilson rela tion approximately have 



been considered. One, the so-called fixed point action (jHasenfratz , 



1998 



, appro ximates the fixed 



Hasenfratz et al, 



point o f a renormalization group transformation by truncating to a small range. 
(|1998|) have shown that (unt runcated) fixed po int fermion actions satisfy the Ginsparg-Wilson re- 
lation. The second version (IGattringen. 120011) . directly minimizes deviations from the Ginsparg- 
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Wilson relation by adjusting the parameters in an arbitrary Dirac operator with a finite (small) 
number of terms. 



C. Numerical simulations 



After having chosen a gauge and fermion action one computes expectation values of interesting 
observables, Eq. ©, by numerical Monte Carlo simulations. For this one creates a sequence of 
gauge field configurations {ujy(x)}, i= 1, . . . ,N, distributed with probability distribution 



pm i \x)}) 



■(detM F (f/)) § exp{-5 G (f/)} 



■exp{-S eff (U)} 



(51) 



z(pT " ^ JJ Z(P) 

Here, 8 = rif, the number of flavors, for Wilson and chirally invariant fermions, and 8 = riff A for 
(rooted) staggered fermions, 5 and now 



S eff (U) = S G (U) - STrlogM F (£/) . 



(52) 



Expectation values (O) are then computed as an average over the ensemble of gauge field config- 
urations, 

(0) = i£o», (53) 



i=l 



where = 0{ujf > ) is the observable evaluated on the gauge field configuration i 

For pure gauge simulations, when no fermions are present, or in the quenched approxima- 
tion, where the fermion determinant is set to one (detMf =1), the action is local (in the gauge 
fields) and the sequence of con figurations can be gene rated with a local updati ng algorithm , 



such as the Metropolis a 



Kennedy and Pendletonl 



gorit 



1985) 



lm ([Metropolis et all 1 1953b or a heatbath algorithm (|Creuta . 



1980; 



With the fermion determinant present, all gauge fields are couple d and the local updating 



rithm s become impractical. Molecular dynamics based algorithms ([Callaway and Rahman . 



algo 



1982. 



1983b have become the standards for simulations with dynamical fermions. For a scalar lattice 



5 The sketch here is somewhat schematic: each fermion with a different mass would get its own determinant factor. 
Furthermore, Mp should be Hermitian and positive semi-definite. For Wilson fermions one therefore takes Mp = 
D' w Dy/ and uses 8 = «//2, while for staggered fermions one takes Mp = [D* KS DKs]ee where the subscript "ee" refers 
to the matrix restricted to the even sublattice. This is possible, since D* ks Dks block-diagonalizes to even and odd 
sublattices. Restricting to only one sublattice removes the doubling introduced by the "squaring." 
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field theory with action S(§ x ) one introduces a fictitious momentum p x on each lattice site, and 
considers the Hamiltonian 



(54) 



(55) 



This Hamiltonian defines a classical evolution in a fictitious time x by, 

; . dS 

where the dot denotes the derivative with respect to x. Given some initial values (p x (0) , <j>-v(0) ) 
these equations of motion define a trajectory (p x (z),§ x (i)) through phase space. The classical 
partition function corresponding to the set of all such trajectories is 



Z = JH[dp x d^ x ]cxp{-H(p^)} = ^ JHd$ x exp{-Sm 



(56) 



where in the second step the quadratic integration over the p x has been carried out, and 9^ is an 
unimportant normalization factor. The integration of Hamilton's equations, Eq. (1551) . conserves 
the Hamiltonian, Eq. (1541) . up to numerical errors. To get the correct distribution corresponding 
to the canonical partition function (1561) , the ficti tious momenta are "refreshed" periodically by 



1985 



1986). This algorithm 



replacement with new Gaussian random numbers (IDuane and KogutL 
goes under the name of Hybrid Molecular Dynamics (HMD). 

Relying on the ergodicity hypothesis, the expectation value of observables can then be com- 
puted by averaging over many MD trajectories 



i r 

{0) = f 



JxO(4»(x)) . 



(57) 



Integration of the equations of motion, Eq. (1551) , is done numerically by introducing a finite 
step size Ax and using a volume-preserving integration algorithm, such as leapfrog. Due to the 
finite step size, the Hamiltonian is not exactly conserved during the MD evolution, leading to 
finite step size errors in observables, including the Hamiltonian itself, of O ((Ax) 2 ) for the leapfrog 
integration algorithm. These step size errors can be eliminated — the algorithm made exact — 
by combin ing the refreshed M D evolution with a Metropolis accept/reject step at the end of each 



trajectory (IDuane et al. 



19871) . resulting in the so-called Hybrid Monte Carlo (HMC) algorithm. 



For a lattice gauge theory the equations of motion have to be set up such that the gauge fields 
remain group elements. This is ensured by writing 



U fl (x)=iH fl (x)U fl (x) 



(58) 
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with Hfj(x) — J\ n t a hj ,{x] a traceless Hermitian matrix and t a the SU(AO generators, see, e.g., 



([Gottlieb et al 



19871) . The MD Hamiltonian is given by 

H(H^x)Mx)) = £ilH£(x) + S eff (U fl (x)) 



(59) 



X,fJ 



The equation of motion for Hjx) is then, somewhat schematically, 

dS eff iU) 



dU„(x) 



(60) 



TH 



where "TH" denotes the traceless Hermitian part. The term on the right-hand side of (1601) is usually 
referred to as the force term. With S e /f of Eq. (l52l) we have 



dS ef f(U) _ dS G (U) 



8Tr 



d£/„(x) 



(61) 



To evaluate (|6TT) we need to know all matrix elements of M F l (U), a dense matrix, even though the 
fermion matrix Mp(U) is sparse. This would be prohibitively expensive. Instead, one estimates 
the inverse stochastically. Let R be a Gaussian random field such that 



R* A (x)R B (y)=b AB b xy , 

where A, B denote color indices, and for Wilson-type fermions also Dirac indices. Then, 

dMf{U) . 



Tr 



-M F \U) 



^-mM^R 



dU^(x 



(62) 



(63) 



and for each random vector R only a single inversion, Mp l (U)R is needed. Typically, for each 
time step in the MD evolution one uses just one Gaussian ran dom vector, and henc e one inversion. 
This algorithm goes under the name of "HMD R-algorithm" (|Gottlieb et a/.lll987f) . 



Instead of doing molecular dynamics starting with S e /f of Eq. (|52l) one can first represent the 
fermion determinant by an integral over bosonic fields, called pseudofermions 



detM F (U) = / ]][d<Z>\x)d<I>(x)]exp{-&M F -\u)<Z>} . 



(64) 



HMD using (1641) . referred to as the 4>-algorithm (|Gottlieb et all |1987|) . consists in creating, to- 
gether with the momenta refreshments, a <3>-field distributed according to Eq. (|64|) 6 and then inte- 
grating the molecular dynamics equations for the effective action 



S eff (U,&) = S G (U) +4> t M F 1 ([/)4> 



(65) 



For Mf = D ' D this can be achieved by creating random Gaussian variables R and then setting <t> = D^R. 
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with the 4>-field fixed. Now the force term becomes 



(66) 



This again requires one inversion, l (U)<i>, in each step of the MD evolution. One major benefit 
of the 4>-algorithm formulation is that an accept/reject Metropolis step is easily implemented at 
the end of each trajectory resulting in an exact HMC algorithm. 

The representation of the fermion determinant by an integral over pseudofermion fields, 
Eq. (1641 ), can formally be extended to fractional powers 8 = n.f/4, as needed for rooted staggered 
fermions, and 8 = «//2, as needed for odd number of flavors for Wilson fermions, 



(detM F (t/)) 5 = f Yl[d&(x)d&(x)]exp{-&Mp 8 (U)&} . 



(67) 



The problem then is, how to deal with Mp . In the HMD R-algorithm this is handled by weighting 
the fermionic contribution to the force by a factor of 8 and evaluating M l R at a point in the 
integration time chosen so that the er rors in observables re main order £ 2 , where £ is the step size in 
the molecular dynamics integration (IGottlieb et all 1 19871) . Clark and Kennedy recently proposed 



using a rational function approximation rewritten as a sum over poles (|Clark and KennedyL 
20051) . 



200. 



4, 



M F 8 (U)^r(M F (U)) = a + 1 £ 



ak 



(68) 



c= [M F (U)+b k 

with suitable constants and b^. A ^-algorithm can then easily be constructed, resulting in 
the so-called rational hybrid molecular dynamics (RHMD) algorithm, or, with inclusion of the 
Metropolis accept/reject step to eliminate errors from nonzero £, the rational hybrid Monte Carlo 
(RHMC) algorithm. Elimination of the noisy estimator yields smaller errors than in the HMD 
R-algorithm at a given integration step size. 

Several improvements of the HMD-type algorithms over the last several years have made them 
substantially more efficient. Th ese improvements include "multiple time step integration schemes" 



(|Sexton and Weingarten , 



ofermion fields ( HasenbuschL 



19921). preconditioning of the fermion d eterminant by multiple pseud 



2001 : 



Hasenbusch and JansenL 



20031). and replacing the leapfrog in 



tegration scheme with more sophisticated "Omelyan integrator s" (|Omelyan et all 



Sexton and Weingartenl 



1992 : 



2002a. 



Takaishi and de Forcrand . 



2006) 



b. 



2003 
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D. Asqtad improved staggered fermions 



Staggered fermions, with only one component per lattice site, and the massless limit protected 
by a remnant even/odd U(l) e xU(l) chiral symmetry, are numerically very fast to simulate. One 
of the major drawbacks is the violation of taste symmetry. At a lattice spacing a of order 0. 1 fm, 
which until recently was typical of numerical simulations, the smallest pion taste splitting Eq. (|39l ) 
for standard staggered fermions is of order A(mj,) = a 2 dp ~ (300 MeV) 2 , i.e., more than twice the 
physical pion mass. Even when the lattice spacing is reduced to about 0.05 fm this smallest splitting 
is still the size of the physical pion mass. It is therefore important to reduce taste violations. 
Since the different taste components live on neighboring lattice sites and in momentum space have 
momentum components that differ by n/a, emission or absorption of gluons with (transverse) 
momentum components close to n/a can change the taste of a quark. Exchange of such ultraviolet 
gluons thus leads to taste violations. 



Su ppressing the coupling to such UV gluons thus should reduce the taste violations (IBlurr 



19971 : 



Lagae and Sinclaii . 



1999 



Lepage! . 



1998 



Orginos and Toussainti 



1999; 



Orginos et al. 



et all 



199 



This can be achieved by replacing the link field in the covariant difference operator V A , (see 
Eq. (fT3l ) by a smeared link built from 3-link staples ("fat3") 



U M (x) -> Uf 3 (x) =T f %(x) = U„(x)+®a 2 £ Ajc/^x) 



(69) 



where the superscript I indicates that the Laplacian acts on a link field, 

Aj[/ jU (x) = (u v (x)U^x + av)U^(x + afj)+U^(x-av)U^x-av)U v (x-av + afj)-2U M (x)) . 



(70) 



In momentum space, expanding to first order in g, Eq. (l69l leads to 



Afi(p) — > Afj(p) + (0 £ {2A A ,(» [cos(apv) - 1] +4sin(a/? /U /2) sm(ap v /2)A v (p)} . (71) 

Choosing co = 1/4 eliminates the coupling to gluons A^(p) with a single momentum component 
Pm = n/a. Adding a 5-link staple ("fat5") 

U,(x)^Uf 5 (x) = !ff%(x) = Uf 3 (x) + ^ £ A £ p A%(x), 



(72) 



eliminates the coupling to gluons with two momentum components p v = %ja and adding a 7-link 
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staple ("fat7") 



(73) 



eliminates the coupling to gluons with all three transverse momentum components p v = %/a. 
For smooth gauge fields, with p m 0, the Laplacian, Eq. (1701) . becomes 



A^[/ /J (x)=aZ) v F V/J + --. 



(74) 



where ■ ■ ■ represent higher order terms in a. The change in Eq. (1691) thus produces a change 
~ a 2 D v F v ^ to the gauge field A^. This is a new (a 2 ) lattice artifact, and wil l occur when u sing 



1999) 



fat3, fat5 or fat7 links. It, in turn, can be canceled by a "straight 5-link staple" (ILepagd . 

A^U^x) = ^(u v (x)U v (x + av)U^x + lav)Uj(x + av + aju)Uj(x + aft 

+ t/J (x - cti) (x - lav) U M (x - lav) U v (x - lav + ap) U v (x-av + aft) - lU M (x] 



= aD v F v ^ + 



(75) 



referred to as the "Lepage-term." In momentum space, expanding to first order in g, this becomes 
— {A^p) [cos(2ap v ) - 1] + 2 sin(a^/2) [sin(ap v /2) + sin(3op v /2)] A v (p) } , (76) 

and thus does not affect the coupling to gluons with momentum components at the corners of the 
Brillouin zone. Therefore, replacing 

,2 



U,(x) - Uf 7L (x) = ffn-U^x) = Uf\x) - U - £ Al%(x 



(77) 



eliminates, at tree level, the coupling to gluons with any of the transverse momentum components 
p v = n /a without introducing new lattice artifacts. 



1989) to 



Finally, for a complete O (a 2 ) improvement we include a so-called "Naik-term" (|Naikl . 
improve the free propagator, and hence the free dispersion relation. To keep the structure of the 
couplings to the different tastes unchanged, this involves adding a 3-hop term, 



v„x(*) - v / ,xW--(v,) 3 xW 



(78) 



—U^(x — ap) (x — lafi) (x — 3ap) %(x — lap) j . 
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FIG. 3 Illustration of taste violations for staggered fermion actions with various link fattenings. The valence 
quark masses were adjusted to give the same m n5 /m p = 0.55 for all fermion actions. The results are for 
quenched gauge field configurations with a Symanzik improved gauge action using p = 7.30. The staggered 
fermion actions are standard, or one-link (OL), "fat3+Naik" (OFN) , "fat5" and "asqtad". The pions are 
labeled by the taste structure, with the taste singlet the hea viest, and the taste ps eudoscalar the pseudo- 



Goldstone boson, the lightest. For more comparisons see (|Orginos et al. 
In the free inverse propagator this changes 



2000). 



-sin(ajc jU ) — > -sinfapp) 



1 2 
l + -sin (ap M ) 



Pv + 0(a 4 ) 



(79) 



The fermion action with only the improvement in Eq. (1791) is referred to as the "Naik action". This 
is also the free (noninteracting) limit of the asq and asqtad fermion actions, defined next. 

We now have all the ingredients for an improved staggered fermion action, called the "asq" 
action (0(a 2 ) improved action): use the covariant derivative with the Naik term, Eq. (1791 , and in 
the one-link term replace the gauge links [/„ by the fat7 links with Lepage term llf 11 of Eq. (TTTT ). 
Replacing the various coefficients in the asq action by tadpole improved coefficients finally gives 
the "asqtad" fermion action. The reduction of taste violations for pions with increasing amount of 
link fattening is illustrated in Fig. |3] 

The Naik term, Eq. (1791) . reduces the lattice artifacts in the pressure for free fermions, and thus 
in the very high temperature limit of QCD as illustrated in Fig. HI left panel, and in the 'sp eed 



Bernard et al. 



1998 ). In 



of light' determined from the pion dispersion relation, right panel, from 
Fig, ffl left panel, "p4" fermions are another variant of improved staggered fermions (|Heller et all 
1999b designed to improve the dispersion relation and high temperature behavior. The speed of 
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FIG. 4 The pres s ure (l eft) per fermion degree of freedom for free Kogut-Susskind, Naik, Wilson and 



p4" (IHeller et al. 



199% fermio ns as a function of N t = l/(aT). The continuum value is shown as th e 



horizontal solid line. Figure from 



Bernard et al. 



(120051) : an earlier version appeared in 



Bernard et al. 



(1998). 



The 'speed of light squared', (rig ht), calculated from the pion dispersion relation, for Naik and K-S pions. 



Figure from 



Bernard et al. 



1998). 



light, shown in the right panel, is determined from pion energies E n (p) for various momenta as 

, El{p)-El{Q) 



(80) 



The O (a 2 ) improvement of the asqtad action gives a sta ggered fermion formu lation with good 



scaling properties, as shown in Fig. [5] for a quenched study (|Bernard et al. 



2000a). 



E. Highly improved staggered fermions 

The largest contribution to the 0(a 2 ) error in the asqtad action originates from the taste- 
exchange interactions. This error can be completely eliminated at one-loop level by adding four- 
quark interactions (which are hard to implement in dynamical simulations) or greatly reduced by 
additional smearings. Multiple smearings, for instance 

U,(x) ^X^x) = ? flL f flL U,(x) (81) 

are found to further reduce mass splittings between pions of different taste. However, they increase 
the number of products of links in the sum for X^x) links and effectively enhance the contribution 
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FIG. 5 Rho masses (left) and nucleon masses (right) in units of r\ « 0.32 fm, in a slight update from 



Bernard et al 



d2000al) . Octagons are unimproved staggered fermions with Wilson gauge action, diamonds 
are unimproved staggered fermions with Symanzik improved gauge action, crosses are Naik fermions and 
squares are asqtad fermions, both with Symanzik improved gauge a ction. For compariso n we also show 



tadpole clover improved Wilson fermions with Wilson gauge a ction (IBowler et al. 



and with Symanzik improved gauge action (iCollins et al. 



2000T) (fancy squares) 



1997!) (fancy diamonds). 



of two-gluon vertices on quark lines (see [Follana et al.\ (120071) for a more detailed discussion). 
Thus, an operation that bounds smeared links needs to be introduced: 



(82) 



where 11 is an operation that projects smeared links onto the U(3) or SU(3) group. Cancellation 
of the 0(a 2 ) artifacts introduced by fat7 smearing with the Lepage term can be achieved on the 
outermost level of smearing, and Eq. (|82l ) can be simplified: 



U,(x) = J^UJ^U^x) = 7 HISQ U,(x) . 

Introducing smeared and reunitarized links that arise after each operation in Eq. (|83l) 



(83) 



V,(x) = 7 fl U,{x) , 
W,(x) = <aV tl {x) = U!F fl U lx {x) , 
X M (x) = f flL W^x) = 9 ms Qu ii {x) 



(84) 
(85) 
(86) 
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we can write the covariant derivative that replaces the naive one: 



Wxto - V,(x)[I]xW--(l+e)(V,) 3 [lV]x(x) . 



(87) 



Equation (187]) is a recently propo sed "Highly improved staggered quark", or "HISQ", discretiza- 
tion scheme (|Follana et all |2007[ ). In square brackets we indicate which links are used as gauge 
transporters in the derivatives. The second term is the Naik term evaluated using the reunitarized 
links V^u(jc). Its coefficient includes a correction £ introduced to compensate for the order (am) 4 
and a s (am) 2 errors. This correction is negligible for light quarks, but may be relevant for charm 
physics if a level of accuracy better than 5-10% is de sired. The correctio n £ can be either tuned 



20071) . 



nonperturbatively or calculated in perturbation theory (IFollana et all 

The HISQ action suppresses the taste-exchange interactions by a factor of about 2.5 to 3 com- 
pared to the asqtad action, which makes it a very good candidate for the next generation of simula- 
tions with 2+1 or 2+1+1 flavors of dynamical quarks, where in the latter case the last quark is the 
charm quark. We discuss preliminary studies of the HISQ action in more detail in Sec.|Xj 



HI. STAGGERED CHIRAL PERTURBATION THEORY AND "ROOTING" 



A. Chiral effective theory for staggered quarks 



Because simulation costs increase with decreasing quark mass, most QCD simulations are done 
with the masses of the two lightest quarks (up and down) larger than their physical values. The 
results, therefore, have to be extrapolated to the physical light quark masses. This is done using 
chiral perturbation theory, the effective field theory that describes the light quark limit of QCD 



(IGasser and Leutwylerl 



1984, 



1985 



Weinberg , 



1979L 



Even with the asqtad improvement of staggered fermions, taste- symmetry violations are not 
negligible in current simulations. It is therefore important to include the effects of discretization 
errors in the chiral perturbation theory forms one uses to extrapolate lattice data to physical light 
quark masses and to infinite volume; in other words, one needs to use "staggered chiral perturba- 



tion theory" (SXPT). Ind eed, it is not 



possib 



continuum chiral forms (lAubin et al. 



e to fit the mass dependence of the staggered data to 
2004bJ). Once the discretization effects are included explic- 
itly by making SZPT fits, one can gain good control of the errors from the continuum extrapolation. 
Furthermore, the effects of taking the fourth root of the staggered determinant can be included in 
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SXPT. The resulting "rooted staggered chiral perturbation theory" (rSXPT) allows us to understand 
the nonlocal and nonunitary consequences of rooting on the lattice and to test that these sicknesses 
go to zero as a — » 0. 



Lee and Sharpel (| 1999b first developed SXPT for a single staggered flavor (a single staggered 



field) a t 0(a 2 ); this was generalized to an arbitrary number of flavors by 



Aubin and Bernard 



2003allb|) . Here, we outline the theory with Nf flavor s to th is order; for the next order we re 



fer the reader to the literature (ISharpe and Van de Watei , 



2005). 



To derive SXPT, one s tarts b y determining, to the desired order in a 2 , the Symanzik effective 
theory (SET) (ISymanzikl . Il983|) for staggered quarks. The SET is an effective theory for physical 
momenta p small compared with the cutoff (p <C 1 /a); it parametrizes discretization effects by 
adding higher-dimensional operators to continuum QCD. In particular, taste violations appear to 
(a 2 ) in the SET as four-quark (dimension six) operators. These operators arise from the exchange 
of gluons with net momenta ~7t/a between two quark lines. Such gluons can change taste, spin, 
and color, but not flavor. Therefore, the operators generated have the form 

Oss'u' = qiiis ® kt)qi qjiis' ® ^)qj , (88) 

where i and j are flavor indices, spin and taste matrices have the notation of Eq. (I3TT) . and color 
indices are omitted because they play no role in what follows. The SU (Nf) vector flavor symmetry 
guarantees that O ss 'tt' is a flavor singlet, which means that i, j are (implicitly) summed over their 
Nf values in Eq. (f88l) . 

The possible choices of the spin and taste matrices in Eq. (1881) are constrained by the staggered 
symmetries. First of all, we can use the separate U(l) e for each flavor. This forces each of the 
bilinears making up O ss i tt i, for example qi(y s <S> ^)qu t0 be U(l) e invariant by itself for each i. 
From Eq. (1371) . we then have that {75 £g) ^5, y s ® = 0, which gives twelve choices for y s and \ t : 
One of them must be a scalar, tensor, or pseudoscalar (S, T or P) and the other must be a vector 
or axial vector (V or A). For example, we might have A <8> T, that is, y s <S) £,t — Y^5 ® ^vX' w i m m e 
notation = 7^75 (and similarly for tastes), and = j^vj^x] ( an d similarly for spins). Such 
operators are called "odd" because, in the original one-component form of Eq. (124)) . the fields % 
and % are separated by an odd number of links (1 or 3) within an elementary hypercube. This is 
easily seen from the equivalence given in Eq. (|3T| ). 

Shift symmetry gives the next constraint. As mentioned following Eq. (l36l) . shift symmetries are 
a combination of discrete taste symmetries and translations. In the SET, however, where external 
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momenta are always small compared with the cutoff, it is possible to redefine th e fields q(y) to 
make t he action invariant under arbitrary translations, like in any continuum theory (|Bernard et all 
2008al) . The shifts then have the form: 



q{y) - (I®^)q(y) 5 fo) - foW® ^) ■ (89) 

Thus, for each of the sixteen possibilities for £, t , the bilinear qi{y s ®^,t)qi undergoes a unique set of 
sign changes under shifts in the four directions //. Since the only bilinears that are invariant under 
all shifts are those with ^ t = I, this immediately shows why taste symmetry cannot be broken by 
bilinear operators. Moreover, it forces \ t — \ t < in the four-quark operators of the SET, Eq. (1881) . 

We now consider the implications of rotations and parity. Rotational symmetry requires that 
Lorentz (Euclidean) indices be repeated and summed over, but since the lattice action is invari- 
ant only under 90° rotations, an index can be repeated any even number of times before sum- 
ming, not just twice. Further, with staggered quarks, the lattice r otational symmetry transforms 



the taste indices together with the space-time (and spin) indices (|van den Doel and SmitL 



Golterman and S mit. 



1983 



1984|) . Since, ^ f = "t, t i, the spin indices on y s > must be the same as those on y s . 
Parity then forces y s and y s i to be identical; combinations such as y s = y v , y s i = y V 5 are forbidden. 
There are now only two choices: either the spin indices and taste indices are sep arately summed 



over, o r there are some indices that are common to both the spin and taste matrices. Lee and Sharp 



dl999J) called the former class of operators "type A," and the latter, "type B." 

Because there are twelve choices for an odd bilinear, there are a total of twelve type- A operators. 
An example is 

0[ V xp] = a z q'i(y»®£>5)qiqj(y t i®%5)qj , (90) 

with the repeated index /j summed over. The fields here have standard continuum dimensions, 
so we write explicit factors of a to give the operator dimension four. Note that type-A operators 
are invariant over the full Euclidean space-time rotation group, SO(4), as well as a corresponding 
SO(4) of taste, a subset of the complete SU(4)y of taste that appears in the continuum limit. 

In order to have a sufficient number of indices to construct a type-B operator, either y s = y s J 
or \ t = \ t i must be a tensor (T); the other set is then either V or A. Thus there are four type-B 
operators. An example is 

°[v ll xT fl ]=a 2 [4KY^®^vk;4y(Y^®^v^ , (91) 
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where the second term ensures that the operator has no separate spin- or taste-singlet piece. Since 
the index /j is repeated four times, one sees explicitly from Eq. (|9T| ) that type-B operators are 
invariant only under joint 90° rotations of spin and taste. 

The SET to O (a 2 ) for Nf flavors of (unrooted) staggered fermions is then simply the continuum 
QCD Lagrangian for ANf species together with the above type- A and type-B operators. 7 Given this 
SET, the (a 2 ) chiral Lagrangian is constructed by finding — with a "spurion" analysis, outlined 
below — the chiral operators that break the full SU(4V/)lxSU(4/V/)a symmetry in the same way 
as the four-quark operators in the SET do. However, the symmetry is also broken by the quark 
mass terms in the SET. In order to arrive at a consistent expansion scheme (a consistent power 
counting) for the chiral theory, we must first decide how the breaking by a 2 terms compares with 
the breaking by mass terms. 

The standard power counting, which we follow here, takes a 2 ~ m, where m is a generic quark 
mass. More precisely, we assume that (see Eq. (1391 ) 

a 2 8 ~ ml p = 2Bm , (92) 

where a 2 d is a typical pion taste- splitting. The taste splittings and squared Goldstone pion masses 
are indeed comparable in current MILC ensembles. Goldstone pion masses range from about 
240 MeV to 600 MeV; while, on the "coarse" (a ~ 0. 12 fm) ensembles, the average taste splitting is 
about (320 MeV) 2 . This splitting drops to about (210 MeV) 2 on the "fine" (a « 0.09 fm) ensembles 
and to about (125 MeV) 2 on the "superfine" (a ~ 0.06 fm) ensembles. It is clear that Eq. (|92|) is 
appropriate in the range of lattice spacings and masses we are working on. However, for future 
analysis of data that include still finer lattices and omit the coarse and possibly the fine ensembles, 
it might be reasonable to use a power counting where a 2 is taken to be smaller than m. 

To derive the leading order (LO) chiral Lagrangian, we start with the Lagrangian in the contin- 
uum limit, i.e., in the absence of taste-breaking operators. In Euclidean space, we have 

2 2 

L cont = L^a^^t) _ l - B f 2 Tx{M E + MI?) + g(7>(d>)) 2 , (93) 

where the meson field <I>, E = exp (/<!>//), and the quark mass matrix M are ANf x ANf matrices, 
and / is the pion decay constant at LO. The field E transforms under SU(4iV/)z. x SU(4ivy)/? as 

7 There are additional O (a 2 ) terms in the SET, for example from the gluon sector, that we ignore here for simplicity. 
Such terms are taste invariant, and at leading order only produce "generic" effects in the chiral Lagrangian: O (a 2 ) 
changes in the physical low energy constants. 
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E -> LE7? t . The field 4> is given by: 



4> 



U Tl+ K+ 
71 D K° 
K K° S 



(94) 



\ ::: 'J 

where each entry is a 4 x 4 matrix in taste space, with, for example, 7t + = X)i=i ^a- The 16 Her- 
mitian taste generators T a are T a = {^s,i^«5, z^/vG" > v),^,7}. Since the normal staggered mass 
term is taste invariant (see Eq. (1321)). the mass matrix has the form M = diag(m J( /, m^I, m s I, ■■ ■). 

The quantity mo in Eq. (|93T ) is the anomaly contribution to the mass of the taste- and flavor- 
singlet meson, the r\' oc Tr(4>) . As usual, the T)' decouples in the limit mn — > °°. However, o ne 
may postpone taking the limit and keep the T|' as a dynamical field ( 



Sharpe and Shoresh . 



200JJ) in 



order to avoid putting conditions on the diagonal elements of 4>. These diagonal fields, £/,/),. 



are, then simp 



dSharpd . 



1990 



the u u, dd bound states, which makes it easy to perform a "quark flow" analysis 



19921) by following the flow of flavor indices through diagrams. 
Since a typical pion four-momentum p obeys p 2 ~ ~ 25m, both the kinetic energy term 
and the mass term in Eq. (l93l) are o(m). By our power counting scheme, Eq. (l92l) . we need 
to add (a 2 ) chiral operators to complete the LO Lagrangian. These are induced by the (a 2 ) 



operators in the SET. We start with the type-A operator 0[ VxP j of Eq. (l90l) . Using q\ = qf + qf, 

R L R Zj 

with q i ' = [(1 ± Y5) /2]qi, and similarly for qi with q t ' = #,-[(1 =FY5) /2], we have 

[VxP] =a 2 [qf {y^ 5 )qf + qf{y^5)qf} 2 = [q* {j, ® F R ) q R + q L (y, ® F L ) q L } 2 , (95) 

where flavor indices are implicit in the last expression. The spurions Fr and Fl will eventually take 
the values 



F R = a ^ } = a ^ 5 ® / flavor , F L = a V ' = a ^5 ® /flavor , (96) 

but for the moment are given spurious SU(4iV/)z,xSU(4iV/)^ transformation properties Z 7 ^ — > 
RFrR^ and — > LF^L} in order to make 0[y x p] "invariant." 

The corresponding o(a 2 ) operators in the chiral Lagrangian are then invariants constructed 
only from E, E^, and quadratic factors in Fr and/or Fl. We cannot use derivatives or factors of the 
mass matrix M because such terms would be higher order. It turns out that there is only one such 
operator: 

. (MA . (MA j. 

(97) 



C 1 Tr(F L EF^E t ) = C^Tv^I^f^ 
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where C\ is a constant that can be determined in principle by fits to staggered lattice data. 

The eleven other type-A operators can be treated in the same way, though of course different 
operators will have different spurions with different transformation properties. Some of the type-A 
operators give more than one chiral operator, but, because of repeats, a total of only eight chiral 
operators are generated. 

The type-B operators couple space-time and taste indices, and are invariant only under 90° rota- 
tions. Their chiral representatives must t herefore have derivatives to carry the space-time indices; 



an example is Tr^E^^E^E^) Jsharpe and Van de Water . 



2005). Because of the deriva- 



tives, however, these operators a re higher order a nd do not appear in the LO chiral Lagrangian. 
This was an important insight of iLee and Sharpd (I1999|) . It means that at LO the physics has the 
"accidental" SO(4) taste symmetry of the type-A operators. 
We can now write down the complete LO chiral Lagrangian: 
a 1 ™2 



L = Lj^IB^) - l -BfTv{M E + M tf) + ^ (Tr(d>) ) 2 + a V 



(98) 



- V 



where the taste- violating potential V is given by 

C lT r(^E^Et) + |[Tr(^^E) +h,.) 
+^[Tr&f^h) +h.c] + f Tr^^) 

_l_ U ^ I _L 



2 

Civ 



4 

Csv 



)+h.c.] + ^[Tr(^ 5 f) m 



[Tr(^ ) E)Tr(^" /; ^)] + ^ [Tr(^E)Tr(^ 



:(*/Vt> 



2 L \-3V — /— V-3V — /J ' 2 L--\-3v5 ^J ±± \^>5V 



(99) 



with implicit sums over repeated indices. 

Expanding Eq. (|98l ) to quadratic order in the meson field 4>, we find, as expected, that pions 
with nonsinglet flavor fall into SO(4) taste multiplets, labeled by P, A, T, V, S. We show numerical 
evidence for this in Sec. IIII.CI The splittings S t of Eq. (l39l) . with t = P, A, T, V, S, are given 
in terms of C\, C3, C4 and C(,. The presence of two traces in the terms multiplied by C2V, C2A, 
C^y, and Csa means that they cann ot contribute at this order to the masses of (flavor-)charged 



mesons. 



Aubin and Bernard! (|2003a|) showed that such terms do generate "taste hairpins," which 



mix the flavor-neutral mesons however, of taste V (and separately, taste A). In other words, there 



a 2 8> 



a 2 8' 



are terms of form —^-(Up + D M + S M H ) and — ^(L^s + + H ) in the expansion 

of Eq. (|98l) , where 8' v and d' A are functions of C2V, Cz4> C$y, and C^a- These terms have been 
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indirectly observed (|Aubin et all l2004bl) in fits to charged pion masses and decay constants to 
one-loop expressions derived from Eq. <l98l) . Because of the practical difficulties in simulating 
disconnected diagrams, taste-hairpins have not yet been studied directly in two-point functions of 
neutral mesons. 



So far, the entir e discussion of SXP ' 



Bernard! (120021) and 



has b een in the context of unrooted staggered quarks. 



Aubin and Bernard! d2003a|) proposed that rooting could be taken into account 



by using quark flow to determine the presence of closed sea-quark loops in an SXPT diagram, and 
then multiplying the diagram by a factor of 1/4 for each such loop. This is a natural assumption , 
becau se it is exactly what happens in weak coupling perturbation theory (|Bernard and Goltermanl 



1994). In the chiral theory, however, the validity of the prescription is not obvious. 



To study in more detail how rooting should be handled in SXPT, it is convenient to replace the 
quark- flow picture with a more systematic way to find a nd adjust the sea-quark loo ps. This is pro- 
vided by a "replica rule," introduced for this problem by lAubin and Bernard! (|2004f) . Since rooting 
is defined as an operation on sea quarks, it is useful first to separate off the valence quarks by 
replacing the original theory with a partially-quenched one: introduce new (valence) quarks along 
with ghost (bosonic) quarks to cancel the valence determinant . The adjustment to the SXPT th eory, 



Eq. (1981) . is the standard one for a partially-quenched theory ((Bernard and GoltermanLll994h : just 
add some additional quark flavors and corresponding bosonic flavors. The masses of the valence 
quarks may be equal to or different from the sea masses. The latter case is clearly unphysical, but 
is useful for getting more information out of a given set of sea-quark configurations. 

We now replicate each sea-quark flavor n r times, where n r is a positive integer, so that there 
are total of n r Np flavors. We then calculate as usual with the replicated (and partially-quenched) 
version of Eq. (l98l) . going to some given order in chiral perturbation theory. The result will be a 
polynomial in n r , where factors of n r arise from summing over the indices in chiral loops. Finally, 
we put n r = 1/4 in the polynomial. We thus take into account the rooting by effectively counting 
each sea-quark flavor as 1/4 of a flavor, which cancels the factor of 4 that arises from the taste 
degree of freedom. The chiral theory obtained by applying this replica rule to SXPT is called 
"rooted staggered chiral perturbation theory" (rSXPT). 

Note that we have done nothing to the valence quarks. Since the number of tastes of the sea 
quarks has been reduced by a factor of 4, it is clear that there is a mismatch, even when the valence 
masses are taken equal to the sea masses. This is still true in the continuum limit, where the issue 
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is particularly transparent. When taste symmetry is exact, rooting removes three of the four tastes 
from the quark sea for each physical flavor, but leaves the valence quarks unaffected. It is therefore 
possible to construct Greens functions, either at the quark or the chiral level, which are unphysical, 
in the sense that the external particles ha ve no counterp art in the intermediate states. Sharpe has 
called this the "valence-rooting problem" dSharpd . l2006bh . The solution is however straightforward 



([Bernard et al. 



2007b 



2008c; 



iharpeL I2006bl) : the physical subspace can be obtained simply by 



choosing all external particles to have a single value of taste (taste 1, say). Using flavor and taste 
symmetries, other Greens function s may also be construc ted that happen to equal these physical 



correlators in the continuum limit (IBernard et al. 



2007bl) . Nevertheless, most Greens functions 



will be unphysical. This is not a cause for concern as long as there is a physical subspace. In fact 
such a situation has nothing, per se, to do with rooting: it will happen in continuum QCD, or in 
any lattice version thereof, if we introduce arbitrary numbers of valence quarks. 

We emphasize that the replica rule tells us to take into account only the explicit factors of n r 
from chiral loops. Putting n r = 1/4 in the polynomial resulting from the SXPT calculation is 
thus a well-defined procedure. We are not concerned with the fact that, if replication is done in the 
fundamental, QCD-level theory, the low energy constants (LECs) such as / and B will be (implicit) 
functions of n r . Such dependence is in general unknown and nonperturbative, and not amenable to 
analytic continuation in n r . Instead, as is always the case in chiral perturbation theory, we treat the 
LECs as free parameters. After setting explicit factors of n r to 1/4 in our calculations, the LECs 
can be determined by fitting the lattice data to the resulting chiral forms. The unknown dependence 
of the LECs on n r is however an obstacle in trying to show, directly from the fundamental theory, 
that rSXPT is the correct chiral theory. This is discussed further in Sec. IIII.CI 



B. Extensions of staggered chiral perturbation theory 



The purely staggered theory discussed thus far is often insufficient for calculations of many 
physical quantities. It would be very difficult, for example, to simulate heavy quarks with the 
asqtad action at currently-available lattice spacings because of the large discretization errors that 
appear when am ~ 1. Thus, the determination of phenomenologically important properties of 
heavy-light m esons and baryon s has u sually been car ried out by adding a heavy v alence quark with 



the Fermilab (El-Khadra et al. 



19971) or NRQCD dThacker and Lepagel 119911) action to asqtad 
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simulations of the sea quarks and light valence quarks. Alternatively, HISQ valen ce quarks have 
been used on the asqtad sea configurations to get precise results for charmed mesons (IFollana et all 



2008b . To the accuracy strived for in current calculations, the effects of heavy sea quarks can be 



neglected; that is, these quarks can be treated in the quenched approximation. 

For several other quantities, the complicated effects of taste- symmetry violation make staggered 
quarks difficult to use. Since these effects often present the greatest o bstacle in the valence sector, a 



very successful compromise, first introduced in 



Renner et al. 



(|2005l) . has been to add domain- wall 



valence quarks on top of the MI 



ing used to st udy scalar mesons (lAubin et al 



C sea-quark ensem bles. Such "mixed- action " simulations are be 



2008 



Renner et al. 



2008b. Bk and related quantities ( Aubin et al 



2009a), nuc 



Detmold et al. 



eon properties (IBratt et al. 



2009; 



Edwards et al. 



20071). hadron spectro scopy ([Edwards et al 



2006a 



2006b; 



Ha gler et al. 



Walker-Loud et al. 



son scattering (IBeane et al. . 



2008a. 



2008c 



dj), and nuclear-physics topics (|Beane et al. 



2007a, 



2007d . 



2008 



2009), me- 



2008b; 



To take full advantage of simulations with heavy valence quarks or mixed actions, it is useful 
to have chiral effective theories that properly include the discretization effects. We briefly discuss 
such theories, starting with the mixed-action case of domain-wall valence quar ks on a staggered 



sea. T he b asic ideas of mixed- acti on chiral perturbation theory were developed in 



Bar et al. 



2004b and 



Golterman et al. 



(2003, 



(|2005r) for the case of Wilson fermions in the sea and chiral fermions in 
the valence sector. By chiral fermions we mean overlap or domain wall quarks, where we assume 
for domain wall quarks that L s is large enough that the residual ma ss is negligible. The extension 
to chiral valence fermions on staggered sea quarks (|Bar et a/.Ll2005|) is then fairly straightforward. 



Features of mixed-action chiral theory that are universa 



the sea-quark action, have been discussed in 



, in the sense that they are independent of 



(120071 . 



2009b). 



Chen et al. 

Because the valence and sea quarks have different actions, a mixed-action theory lacks the 
symmetries that would normally rotate valence into sea quarks (or vice versa) in a standard theory. 
Since we assume that both the valence and sea sectors approach the expected continuum theories 
as a — > 0, these symmetries should be restored in the continuum limit. At the level of the Symanzik 
effective action, the violation of these symmetries first appears at O (a 2 ) in the existence of inde- 
pendent "mixed" four-quark operators: in our case, the product of a domain- wall (valence) bilinear 
and a staggered (sea) bilinear. We know, following the development in Sec. IIII.Al that each bilinear 
must be separately chirally invariant, and that any staggered bilinear must be taste invariant. It is 
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then simple to see that only two mixed four-quark operators are possible: 



A = a 2 \\f a yfl5Va qi{y5ln®I)(li 



(100) 



where \\f a is a domain-wall quark or ghost of valence flavor a, and qi is a staggered quark of sea 
flavor i, and a and i are summed over. As in the pure staggered case, the color indices in these 
operators are irrelevant. 

In addition to the operators in Eq. (1 1001) . there are the full complement of standard, purely stag- 
gered four-quark operators in the sea sector, and standard, purely domain- wall four-quark operators 
involving valence quarks and valence ghosts. In a normal theory, the relative coefficients of cor- 
responding sea-sea, valence-valence, and valence-sea operators would be fixed by the symmetries. 
But in the mixed case, all such operators are independent and must be treated separately. 

In the corresponding chiral effective theory, the purely sea-quark sector is the same as the 
sea-quark sector of a standard staggered theory. Similarly, the purely valence-quark sector is the 
same as the valence-quark sector of a standard domain- wall theory. Mixed valence- sea mesons are 
affected by various operators, including the operator corresponding to Eq. (11001) : 



« 2 CMixTr(x 3 Ex 3 E 



(101) 



where E is the complete chiral field involving both sea and valence (and ghost-valence) quarks, 
and X3 is a diago nal matrix that takes the value +1 in the sea sector and — 1 in the valence sector. 



At LO one finds ( Bar et al . 



2005; 



Chen et all 



2009a) 



111 



TZ,ab 



111 
111 



%,l],t 
2 

n.ia 



B(m a + m b ) 
B(mj + m.j) +a 2 5f 
B(mi + m a ) +a 2 8 M ix 



(102) 



where a,b are domain-wall (valence) flavors, i,j are staggered (sea) flavors, t is the taste of a 
sea-sea meson, as in Eq. (1391), a nd S m^ is a func t ion of Cyiix an d other low energy constants. 



(|2008|) have determined Smix numerically by 



Orginos and Walker-Loudl (|2008|) and lAubin et al. 
measuring the masses of mixed mesons. 

The mixed-action chiral Lagrangian th us developed can be used to calculate one-loop effects in 



pseudoscalar ma sses and decay con stants (IBar et all 



K — K scattering (|Chen et all 



2006). 



20051) . in B K (lAubin et al 



2007b ) and / = 2 
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Next, we consider the case of heavy-meson staggered chiral perturbation theory (HMSXPT), the 
relevant chiral theory for a heavy meson made out of a heavy valence quark and a light staggered 
valence quark, on the background of staggered sea quarks. HMSXPT is designed to parameterize 
the light quark chiral extrapolation and the light quark discretization effects. Discretization errors 
due to the heavy quark are not included; it is assumed that they can be estimated in dependently by 



using heavy-quark e ffective theory (HQE r 
lattice heavy quark ( Kronfeldl 



19921 : iNeubert . 



1994 ) to describe the 



?) (|Isgur and Wisd . 
|200fl[2004j). 

At the level of the SET, the first nontrivial effect of combining the heavy quark with the stag- 
gered theory is again the generation of mixed four-quark operators (a heavy-quark bilinear times 
a light-quark one). As before, such operators do not break taste symmetry. Furthermore, unlike 
the mixed-action case, symmetry between heavy and light quarks is already strongly broken by the 
heavy-quark mass. So the mixed operators have no important effects in this case. 

The power counting for he avy-light mesons in XPT makes the HMSXPT at LO rather sim- 



ple (lAubin and B ernard 



(|Manohar and Wise . 



20061). I n the continuum, the chiral Lagrangian for heavy-light mesons 



20001) starts at 0(k), with k the residual momentum of the heavy quark. The 



light meson momentum p should also be o(k). In our staggered power counting, Eq. ( |92| ), we 
take p 2 ~ ~ a 2 . This means that the LO heavy-light meson terms are lower order than the 
O (a 2 ) discretization errors in the light quark action. The LO heavy-light meson propagator and 
vertices are thus the same as in the continuum, as are the heavy-light currents that enter, e.g., 
in leptonic and semileptonic decays. The light-quark discretization errors in heavy-light me- 
son quantities first appear at one loop (NLO), through the taste violations in the light meson 
propagators in the loop. These correcti ons have been calculated for heavy-light leptonic de- 



cay constants (lAubin and Bernard . 



( Aubin and Bernard 



(ILaiho and Van de Waten. 



20061) . for semileptonic heavy-to-light decays, e.g., B — > %, 



2007), and for semileptonic heavy-to-heavy decays, e.g.,B — > D and B — > D* 



2006). There are also analytic NLO corrections to physical processes, 



coming both from light-quark mass corrections (as in the continuum) and from taste- violating cor- 
rections to the LO Lagrangian and currents. In practice, it is usually easy to guess these analytic 
NLO corrections from s ymmetry arguments, so it is not necessary to use the complicated NLO 



heavy-light Lagrangian (jAubin and Bernard . 



2006 
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C. The issue of rooting 



The extra tastes are eliminated in staggered dynamical simulations by taking the fourth root of 
the fermion determinant — the fourth-root procedure. In the past few years there has been progress 
in understanding and validating this procedure, and we give a brief overview of t his progress here . 



For more detai 



Kronfeld! (120071) and 



ed di scussion, and ful 



GoltermanlfcOOSI) 



lists of references, see recent reviews by ISharpd (I2006b|) . 



The fourth-root procedure would be unproblematic if the action had full SU(4)y taste symmetry, 
which would give a Dirac operator that was block-diagonal in taste space. Indeed, this is the 
"cartoon version" of what we expect in the continuum limit. Assuming taste symmetry is restored, 
the positive fourth root of the positive staggered determinant would then become equal to the 
determinant of a single continuum species. 

However, at nonzero lattice spacing a, taste symmetry is broken and the Dirac operator is not 
block-diagonal (see Eq. (l32l)). From Eq. (|34|) . one has 



\ndet(D KS + m®I) = 4 In det(Z> + m) +lndet{/+ [(D + m) 1 ®/]aA} . 



(103) 



Since (D + m)~ l is nonlocal, we should no t expect the rooted theory to be local for a ^ 0. In fact it 
is possible to prove ( Bernard et all. 2006b ) that the fourth root of the determinant is not equivalent 
to the determinant of any local lattice Dirac operator. 8 The idea of the proof is simple: If there 
were such a local operator, then one could construct a theory with four degenerate quarks, each 
one with that local action. Calling this introduced degree of freedom "taste," one now has a local 
theory with exact SU(4V taste symmetry by construction, and whose determinant is equivalent 
to that of the original staggered theory. This is a contradiction, because the taste symmetry of 
the constructed theory guarantees that it has fifteen pseudo-Goldstone bosons (pions), whereas the 
stagge r ed pions are known to split up into nondegenerate irreducible representations (IGoltermanl 



1986b; 



Lee and Sharpd . 



19991) . Indeed, Fig. [6] shows our lattice measurements of the pion split- 



tings as a function of quark mass (left) and lattice spacing (right). The left plot clearly shows the 
characteristic splitting of the charged pion multiplet into the five nondegenerate submultiplets 
with tastes P, A, T, V, S. This is as predicted at O (a 2 ) in the chiral expansion, as discussed in 



"Equivalent" here means equal up to a factor of the exponential of some local effective action of the gauge field. 
This is e nough to guarantee that the two theories have the same physics at distances much larger than the lattice 



spacing dAdams , 



2005 



Shamir 



2005). 
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FIG. 6 Squ a red charged p ion masses, in units of r\ , as fun ction of quark mass 



Bernard et al. 



(2006e 



Bernard et al. 



(left). Figure from 



(I200lh . The splittings ap- 



200711) . A previous version appeared in 
pear to be independent of the quark mass. The taste splittings as function of a 2 a 2 (right) in a log-log plot, 
showing the expe cted behavior, i ndicate d by the diagonal straight line. A slightly different version of this 



figure appeared in 



Bernard et al. 



(I2007dl) . 



Sec. IIII.Al Further splitting at 



ligher order into a total of eight submultiplets is allowed by the 



1986bl) . but we see little evidence of that at the current level of 



lattice symmetries (IGoltermanL 
statistics. 

The same features of the rooted theory th at imply nonlocality also imply nonunitarity on the 



lattice ([Bernard . 



2006; 



Bernard et al. 



2007allbl : 



Prelovsek . 



2006bh . The issue is particularly sharp 



in the rooted one-flavor theory. The physical one-flavor theory should have no light pseudoscalar 
mesons (pions) but only a heavy r\'. In a rooted theory with exact taste symmetry (e.g., with four 
copies of rooted overlap quarks), this works automatically: the fourth power of the fourth root of 
a (positive) determinant is equal to the determinant itself. Alternatively, one can check directly in 
the rooted four-taste theory that, in physic al correlators, the pion intermediate states cancel and 



only the r\' remains ([Bernard et al 



2007bl) . On the other hand, in the rooted one-flavor staggered 



theory, the pions have different masses at nonzero lattice spacing and cannot cancel, leaving light 
intermediate states with both positive and negative weights. This is a clear violation of unitarity. 

In the continuum limit, we expect that all the pions become degenerate. For the tree-level 
improved asqtad fermions, generic lattice artifacts are of order (a s a 2 ). Taste violations, however, 
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require exchange of at least two UV gluons, since the coupling of a quark to a single gluon with any 
momentum components equal to n/a vanishes. Therefore taste violations with the asqtad action 
should vanish as a? s a 2 as a — > 0. The lattice-spacing dependence of the pion splittings, shown in 
the right-hand plot of Fig. [6} agrees very well with this expectation. Note that since we are looking 
here at flavor- nonsinglet pions, the taste-singlet Tlf also becomes degenerate with the other fifteen 
pions as the continuum limit is approached. 

Thus, the rooted staggered theory is inherently nonlocal and nonunitary at nonzero lattice spac- 
ing, but should become local and unitarity in the continuum limit if taste symmetry is restored. This 
is because, in the limit of exact taste symmetry, rooting of the sea quarks is equivalent to restriction 
to a single taste, which is a local operation. Clearly, the numerical evidence for taste-symmetry 
restoration in the continuum is strong, and accords with the theoretical expectation coming from 
the fact that taste violation is due to an operator with dimension five. How, then, could rooting 
go wrong? The main problem is that the theoretical expectation is based on standard lore of the 
renormalization group (RG) that operators with dimension greater than four are irrelevant in the 
continuum limit. This standard lore for the scaling of operators assumes a local lattice action, 
which does not apply here. The numerical results indicate that the lore is not leading us astray, but 
of course numerical evidence does not constitute a proof. 

There is a further problem in the formal argument we have made so far that rooting is equivalent 
in the continuum limit to restriction to a single taste. The argument seems to require that taste sym- 
metry is restored for the Dirac operator D K $, Eq. (|34|) . itself. In Fig.|6l however, we are only testing 
the restoration of taste symmetry at physical scales, those much larger than the lattice spacing. At 
the scale of the cutoff, there is actually no reason to ex pect that tas t e symmetry is restored. Ind eed, 



2004; 



Follana et all 



2004) find 



direct studies of the eigenvalues of Dks on the lattice (jDurr et all 
only approximate quartets of eigenvalues (indicating approximate taste symmetry) for low-lying 
ei genvalues, those corre sponding to long (physical) distance scales. 



Shamir (2005 



2007|) has set up an RG framework for both unrooted and rooted staggered theo- 



ries, and used it to address the potential problems of rooting. The renormalization group is clearly 
the natural framework to study the scaling of operators, and it also makes possible a more precise 
treatment of the continuum limit. As one blocks Dks to longer distance scales, the eigenvalues at 
the scale of the cutoff are removed, and one may then expect that taste symmetry is truly restored. 
Shamir's RG scheme starts with unrooted staggered quarks, and blocks them on the hypercubic 



47 



lattice by a factor of 2 at each step, integrating out the finer quark fields. The gauge fields are also 
blocked, but the integration over them is postponed until the end, so that the quark action stays 
quadratic at every step. The starting "fine" lattice spacing ay is blocked n times to a final "coarse" 
lattice spacing a c . As n is increased, the coarse spacing is held fixed but small, with a c <C 1 /Age/)- 
The fine lattice spacing thus obeys a/ = 2~ n a c , and the continuum limit is n — > °°, which sends a/ 
to zero. In this unrooted theory, the scaling of A like aj is guaranteed by the standard lore, since 
the action is local. 

The rooted theory cannot be blocked in the same way because rooted quarks are not defined by 
a standard Lagrangian, but by a rule to replace the fermion determinant by its fourth root in the 
path integral. We can, however, apply the rule at every stage of the (unrooted) blocking, obtaining, 
at the n th step, the theory given by 

z KSroot = J dA det \ + ^ ^ ;) ^ (104) 

where D K s, n is the blocked staggered Dirac operator, m n is the (renormalized) mass on the blocked 
lattice, and dR is the full gauge measure (which includes integrals over gauge fields at each level 
of blocking, as well as Jacobian terms coming from integrating out the fermions on the coarse 
lattices). This defines a RG for the rooted theory. However, it is difficult to make progress directly 
from Eq. (11041) . because of the problem of nonlocality. 

Shamir's key insight is that one may define, at each stage of blocking, an intermediate, 
"reweighted theory," which becomes closer and closer to the rooted staggered theory but retains 
locality. Define D n to be the taste- singlet part of Dks.h, and cifA n to be the remainder: 

D n = ^tr ts (D KS , n ) , 
D K s,n = D n ®I+a f A n , (105) 

where tr^ is the trace over taste, and / is the identity in taste space. This parallels Eq. (l34l) . We will 
see below the explicit af in the second term of Eq. (11051) does not mislead us about the scaling of 
a/A. The operator D n is local because Dks is. Further, det(D n +m n ) = det 1//4 ((D„ +m n ) <g)I). The 
(rooted) reweighted theory is then defined by 

z reweighted = j da det(D„ + m n ) , (106) 

Now, since the reweighted theory is QCD-like, albeit with a more complicated gauge integration 
than usual, we expect it to be renormalizable and asymptotically free. The running of the operator 
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a/A„ from ay to a c can then be calculated perturbatively because in this range the lattice spacings 
are all much less than \/Aqcd- Because the theory is local, the standard lore tells us that the 
perturbative running will be a reliable guide to the complete, nonperturbative behavior. Thus we 
expect that the operator norm of a fA n will obey, in an ensemble-average sense, 

ai 



a f 

IM„II<4 



(107) 



where the < sign implies that the scaling is true up to logs. For the same reasons, the mass m n 
should run logarithmically, just as in QCD. From this and Eq. (11051) . we have 



det*{D KS , n + m n ®I) = det(D n +m n )exp (J trln [/+ ((D n + m„) l ®I)a f A„]) 

a f 



det(D n + m„) 1 + 



(108) 



where the quark mass provides a lower bound to the absolute value of the eigenvalues of D n + m„. 
Thus, 

lim Zf Sroot = lim Z* wei s hted . (109) 

In other words, the nonlocal rooted staggered theory coincides with a local, one-taste, theory in the 
continuum limit, as desired. 

Note that Eq. (11081) makes it clear that one must take the continuum (a/ — > 0) limit before 



the chir al (m 



known ([Bernard . 



0) limit for rooting t o work . This is not surpri s ing, s i nce it is already we 



2005 



Bernard et all 



2007b 



Purr and Hoelbling . 



2005 



Smit and VinkL 



1987) 



that the two limits do not commute for all physical quantities, and that taking the chiral limit 
first can give incorrect answers. This is true even for the unrooted staggered theory. As a triv- 
ial example, consider the low energy constant B (see Eq. (l39l ) defined at a given lattice spac- 
ing a by B(a) = m\ t /(2m) for some taste t. Unless t = P, giving the Goldstone pion, one has 
lim^o lim,„^o B(a) = °°; while the desired result is lim,„^o Hm a ^o B(a) = B. 

The reader may worry that the argument thus far presumes too much about how perturbation 
theory works in the reweighted theory. After all, the perturbation theory involves multiple levels 
of gauge integrations, m aking it quite co mplicated. Indeed, no such perturbative calculations have 
been performed to date. IShamirl (|2007l) has pointed out, however, that we may avoid the details 
of perturbation theory in the reweighted theory by leaning a bit more on the standard lore and on 
perturbation theory in the unrooted staggered theory, which is fairly well understood — see lSharpe 
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(I2006bl) and the references therein. One starts by considering the unrooted staggered theory repli- 
cated n r times, where n r is an integer. In this theory the |3 function and the logarithmic anomalous 
dimension of a,fA n will be the standard functions of the total number of fermion species, and a/A n 
will scale as expected as long as n r is not so large that asymptotic freedom is lost. 

Now, a/A n is just the difference between the (replicated) unrooted staggered theory and a (repli- 
cated) unrooted reweighted theory defined by the Dirac operator (D n + m n ) ®7. Since a/A„ gets 
small as n — ► °° in one theory, it must get small in the other theory. Both theories are local, so the 
standard lore says that a/A n scales as expected in perturbation theory in the unrooted reweighted 
theory — however complicated such calculations would actually be in practice. The results of 
perturbation theory to any fixed order are polynomial in n r , with the power of n r just counting the 
number of closed quark loops. So in this perturbation theory, we may put n r = 1/4 to obtain the 
perturbation theory for the rooted reweighted theory, Eq. (11061) . Thus we do not have to calculate 
explicitly in either the unrooted or rooted reweighted theories; we know that atA n will scale to 
zero as expected in perturbation theory. Now the standard lore takes over, as above, for the local, 
rooted reweighted theory, and says afA n will scale to zero as n — » q° even nonperturb atively. 



Bernard et al. 



(12006d) . The results 



A numerical test of the scaling of o/A n was attempted in 
were encouraging but far from conclusive, due to quite large statistical errors. 

Of course, although the above arguments make it plausible that rooting works, they do not 
constitute a rigorous proof. As always in lattice QCD, one relies heavily on the standard lore 
about RG running of irrelevant operators, which is what "guarantees" universality. Further, we are 
unable to do justice here to all the arguments and assumptions involved in the perturbative analysis. 
We have also ignored the nontrivial issues involving the Jacobian obtained by integrating out the 
fermions at each level of blocking. The Jacobian can be written as the exponential of an effective 
action for the gauge fields. The claim is that this effective action is local, bas ically because it comes 



from short -distanc e fluctua t ions of the fermi ons. The reader i s urge d to see 



reviews by lSharpd (I2006b|) . 



Kronfeldl (120071) and 



Shamirl ( 20071) and the 



Goltermanl (12008b for details and discussion. 



We now consider the question of whether rSXPT is the correct chiral theory for rooted staggered 
QCD. This is important first of all because rSXPT allows us to fit lattice data and take the limits 
a^0 and m — > in the correct order and with controlled errors. In addition, the validity of rSXPT, 
coupled with the strong numerical evidence for the restoration of taste symmetry for a — > (see 
Fig. [6]), guarantees that rooted staggered QCD produces the desired results for the pseudoscalar 
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FIG. 7 Relative weights (shown at the right of each line) of two-particle intermediate states in the scalar, 
taste-singlet correlator in the one-flavor case. The two-r)^ state (5 indicates taste singlet) is shown at top; 
while the various two-pion states below are labeled by the pion taste (S, V, T, A, P). The height of each line 
represents, qualitatively, the relative mass of the state. 

meson sector in the continuum limit. This is because rSXPT becomes continuum XPT when taste 
symmetry is restored. 

Before discussing the arguments supporting rSXPT, we note that rSXPT has the main features 
desired for a chiral effective theory of the rooted theory. In particular rSXPT reproduces the nonuni- 
tarity and nonlocality of rooted staggered QCD at nonzero lattice spacing. This comes about be- 
cause rSXPT, like the rooted staggered theory itself, is not an ordinary Lagrangian theory, but a 
Lagrangian theory with a rule. For rSZPT the rule is: calculate in the replicated theory for integer 
n r number of replicas, and then set n r =1/4. Setting n r = 1/4 gives "funny" relative weights for 
different diagrams, which can result ultimately in negative weights for some intermediate states 
in an ostensibly positive correlator. For example, Fig . [7] shows the weights of various two- meson 



intermediate states coming from a rSXPT calculation (|Bernard . 



2006; 



Bernard et all 



2007a) of the 



scalar, taste-singlet correlator in a one-flavor rooted staggered theory. The physical theory should 
only have a two-T)' intermediate state, but here we have various light pion states, with the taste- 
singlet pions 9 having a negative weight. In the continuum limit, however, all the pions become 
degenerate, and they decouple since their weights add to zero. 
The first argument for the validity of rSXPT was given by 



Bernard 



(120061) . The starting point is 



the observation that the case of four degenerate flavors of rooted staggered quarks is particularly 



9 The taste-singlet pion is distinct from the n/ here because it is a flavor nonsinglet arising at the arbitrary, integral n r 
values at which the calculation is done. 
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simple because it is the same as the case of on e flavor of unrooted stag gered quarks. Thus we know 
the chiral theory: it is exactly that obtained by LLee and Sharpd (I1999Q for one unrooted flavor. This 
chiral theory is equivalent to that of rSXPT for four degenerate flavors. The equivalence is manifest 
order by order in the chiral theory: Since the result for any physical quantity is polynomial in the 
number of degenerate flavors, taking Ahr degenerate flavors and then putting hr = 1/4 gives the 
same chiral expansion as a one-flavor theory. 

The case of four nondegenerate flavors may then be treated by expanding around the degener- 
ate limit. The expansion is however somewhat subtle. Once we move away from the degenerate 
limit, nontrivial weighting factors of various diagrams, caused by the fourth root of the deter- 
minant of the sea quarks, come into play. This means that it is impossible to write all needed 
derivatives with respect to the quark masses as derivatives in the one-flavor unrooted theory of Lee 
and Sharpe. The solution is to keep the sea quarks degenerate, but to introduce arbitrary numbers 
of valence quarks. Bernard then shows that it is possible to rewrite all derivatives with respect to 
sea quark masses as sums of various combinations of derivatives with respect to the valence quark 
masses. This approach allows us to remain in the degenerate sea-quark limit, where the chiral the- 
ory is kno wn. It is however necessary to assume that partially-quenched chiral perturbation theory 



(PQXPT) (jBernard and GoltermanL 119941) is valid in the unrooted case. Since the unrooted case is 
local, this is very plausible. Further, there is a significant amount of numerical work that supports 
the validity of PQXPT for local theories, using other fermion discretizations, not just staggered 
quarks. But it should be pointed out that partially-quenched chiral t heories rest on shakier ground 
than the standar d chiral theory fo r QCD, as emphasized recently by lSharpel (I2006ah . For example, 
the argument by IWeinbergl (| 1 9791) for QCD in vokes unitarity, w hich partially-quenched theories do 



not have. On the other hand, the argument by 



Leutwylen ( 



19941) emphasizes cluster decomposition 



instead of unitarity and may be possible to apply to a p artially-quenched Euclidean the ory. Work 



on putting PQXPT on a firmer foundation is in progress ([Bernard and Golterman , 



2009). 



An additional, technical assumption for this approach is that the mass expansion does not en- 
counter any singularities. This is reasonable because the expansion is about a massive theory, and 
one therefore does not expect infrared problems. 

To reach the more interesting case of three light flavors, Bernard raises the mass of one of the 
four quarks (call it the charm quark, with mass m c ) to the cutoff, decoupling it from the theory. This 
requires an additional technical assumption, arising from the fact that there is a range of masses, 
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which begins roughly at m c ~ 2m s (with m s the strange quark mass), where the charm quark has 
decoupled from the chiral theory, but not yet from the QCD-level theory. While the resulting three- 
flavor chiral theory has the same form as that of QCD when a — > 0, the assumption does leave open 
the possible "loophole" that the LECs have different numerical values from those of QCD. 

The above argument takes place entirely within the framework of the chiral theory. It has the 
nice feature that the recovery of the correct QCD chiral expressions, and the vanishing of nonlo- 
cal and nonunitary effects, only requires taste violations to vanish in the continuum limit in the 
unrooted, and hence local, theories with integral n r . The vanishing of these taste violations in the 
rooted chiral theory then follows. On the other hand, because the argument does not connect rSXPT 
to the QCD-level rooted staggered theory, the replica rule ends up emerging rather mysteriously. 
The chain of reasoning also depends on several technical assumptions. 

An argument for the validity of rSXPT directl y from the fundament al rooted staggered theory 



is therefore desirable. It has been developed by iBernard et all (|2008ah by starting from the RG 



framework of Shamir. The basic idea is to generalize the fundamental (lattice-level) theory to one 
in which the dependence on the number of replicas n r is polynomial to any given order in the 
fine lattice spacing cif. Then we can find the chiral theory for each integer n r in a standard way 
(because the theories are local), and apply the replica rule to get the rooted staggered theory at the 
fundamental level and rSXPT at the chiral level. 

For simplicity we focus on a target theory with n s degenerate quarks in the continuum limit. 
Unlike the previous argument, the extension here to quarks with nondegenerate masses is straight- 
forward. Consider Eq. (11041) . the rooted staggered theory at the n th step of blocking, but with n s 
degenerate staggered flavors: 



yKSroot/ 



[n s ) = JdR det'4 (D KS} n + m n ®I) , (110) 
Now generalize this, using the definitions of Eq. (11051) . to 

J ae,t n r[{D n + m n )®I\ 

where t is a convenient interpolating parameter. When t = 1 and n r = n s /4, this reduces to Eq. (II 101) 
because the determinants of the reweighted fields (those involving D n + m or (D n + m)®I only) 
cancel, and the remaining determinant is just that of the rooted staggered theory. When t = 0, on 
the other hand, Eq. (11111) gives a local theory of n s reweighted one-taste quarks. 
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Equation (II 1 II) has an important advantage over Eq. (II 101 ). While the dependence on n s is 
unknown and nonperturbative in both cases, the dependence on n r of Z,f en (n 5 , n r ) is well controlled 
because it vanishes when the taste violations vanish (a/A„ = or t = 0). This makes it possible to 
apply a replica rule on n r at the fundamental QCD level. To see this, we first write 

det nr \(D n +m n ) <g> I+taf AJ , r s i N 1X 

de,M(fi,+ m „)«/] =expKtr ln [ 1 + ((D„ + m „)-^/) to/ A„]) . (112) 

We then expand in powers of the fine lattice spacing ay. These can come from the explicit factor 
a/ in the taste-violating term or from the implicit dependence on a/ of the gluon action and the 
lattice operators D n and A„. The parameter t serves to keep track of the explicit dependence; the 
power of t must be less than or equal to the power of ay to which we expand. From Eq. (II 121) . 
the power of n r must in turn be less than or equal to the power of t. Thus, to any fixed order in 
a f, the dependence of the theory on n r must be polynomial. This means that n r is a valid replica 
parameter of the fundamental theory (again to any fixed order in af). We can in principle find 
the polynomial dependence of any correlation function by calculations for integer values of n r 
only, and then determine the correlation function in the rooted staggered theory by simply setting 
n r = n s /A (and t = 1). 

We now discuss the effective theories, the SET and the chiral theory. For convenience, we can 
work at t = 1 . For n r and n s (positive) integers, Z,f en (n s , n r ) is a local, but partially-quenched, theory 
that can be written directly as a path integral. It is partially quenched because the determinant in 
the denominator needs to be generated as an integral over ghost (bosonic) quarks. Finding the SET 
and the chiral effective theory for such local theories is standard, although the caveats about the 
foundations of PQXPT apply. All that we really need to know is that the effective theories exist for 
any integer n r and n s , and that their dependence on n r is polynomial (because the dependence in 
the underlying theory is polynomial). In the chiral theory we can then set n r = n s /4. At the QCD 
level this just gives the rooted staggered theory for n s flavors. At the chiral level, the reweighted 
parts of the theory again cancel order by order at n r = n s /4, because we have n s flavors of one-taste 
quarks and n r flavors of four-taste ghost quarks, with exact taste symmetry. We are then left with 
exactly the result we would have gotten from rSXPT. 



Th is argument avoids the "loophole" and technical assumptions of the argument of 



Bernard 



(120060 . It also makes clear how the repl ica rule arises from the fundamental theory. On the 



other hand, it inherits the assumptions of 



Shamni (|2007|) . since it is based on that framework. 
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Both arguments rely on the standard PQ5CPT for local theories. This is not surprising since rooted 
staggered QCD inherently shares some features of a partially-quenched theory: Since rooting is 
done only on the sea quarks, there is an excess of valence quarks. As noted earlier, however, this 
"valence-rooting" issue is not a fundamental problem because there is a physical subspace. 

A nice feature of the current argument is that, by coupling rSXPT directly to the RG framework, 
it makes numerical tests of rSXPT into tests of the RG framework, and hence of the validity of 
rooting at the fundamental level. We discuss such tests in Sec. fVIl 



We now turn to the objections raised to rooted st aggered 



2008 allb|). Since th ese objections have been refu ted. (I Adams , 



Golterman, 



quarks 



by. 



Creutz (2006a. 



20081: iBernard et al. 
2008) — see also the reviews i06bT and 



Kronfeld (2007 



2007a.b 



2007 



b. 



20081 



c. 



we give 



only a very brief discussion here. The main point is that most of Creutz's claims apply equally 
well to the proposed continuum limit theory of rooted staggered quarks: a rooted fo ur-taste theory 



Bernard et al. 



with ex act taste symmetry, which is called a "rooted continuum theory" (RCT) by 
(|2008ch . Such a theory provides a tractable framework in which to examine Creutz's claims. 
Because, as emphasized before, det 1 / 4 ((D + m) ®/) = det(D + m), the RCT is clearly equivalent 
to a well-beh aved, on e -taste theory, and gives a counterexample to most of Creutz's objections. 



Alternatively, 



Adams 



(120081) has found counterexamples to Creutz's claims in a simple lattice 
context, namely a version of twisted Wilson quarks. 

While the RCT is equivalent to a one-taste theory, it is not exactly the same in the following 
sense: In the RCT, with its four tastes, one can couple sources to various tastes and generate Green 
functions that have no analogue in the one-taste theory. Such unphysical Greens functions are at 
the basis of many of the "paradoxes" Creutz finds. For example, one can find 't Hooft vertices that 
are singular in the limit m — > 0. Nevertheless these unpleasant effects exist purely in the unphysical 
sector of the RCT; in the physical sector all 't Hooft vertices are well behaved. 

Finally, Creutz has noticed that there is a subtlety involving rooted staggered quarks for negative 
quark mass, and this is in fact true. Independent of the sign of the quark mass, the staggered deter- 
minant is positive, as discussed following Eq. (1281) . The fourth root of the determinant generated 
by the dynamical algorithms, Sec. Ill.Ci is then automatically positive for any sign of m. In other 
words, the rooted staggered theory is actually a function of |m|, not m. This means that rooted 
staggered fermions cannot be used straightforwardly to investigate the effects that are expected 
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(iDashenl . 1 197 11 : IWittenl 1 19801) to occur for an odd number of negative quark masses. 10 A related 
problem occurs when one adds a chemic al potential to the theor y — the determinant becomes 
complex, and the fourth root, ambiguous (IGolterman et all 120061) . Nevertheless, these problems 
have no relevance to the validity of the rooted staggere d theory in the usual ca se of positive quark 



mass and zero chemical potential. For more details, see 



Bernard et al. 



(2007b). 



IV. OVERVIEW OF THE MILC LATTICE ENSEMBLES 



In this program of QCD simulations, ensembles of lattices were generated at several lattice 
spacings and several light quark masses. This allows extrapolations to zero lattice spacing (the 
"continuum extrapolation") and to the physical light quark mass (the "chiral extrapolation"). In 
all ensembles the masses of the up and down quarks are set equal, which has a negligible effect 
(< 1%) on isospin-averaged quantities. The original goals of the program were to simulate with 
three dynamical quark flavors, with a large enough physical volume to make finite size effects 
small, and to vary the quark masses to study the effects of "turning on" the dynamical quarks. 
It later became clear that more lattice spacings were needed to understand the continuum limit. 
Fortunately, computer power was increasing rapidly, which made the simulations with smaller a 
practical. 

Currently, the lattice spacings of the ensembles fall into six sets, with lattice spacings approxi- 
mately 0.18 fm, 0.15 fm, 0.12 fm, 0.09 fm, 0.06 fm and 0.045 fm. In many places these are called 
"extracoarse," "medium coarse," "coarse," "fine," "superfine," and "ultrafine," respectively. The 
0.12 fm lattices were the first to be generated. Over time, as computer power permitted, the lattice 
spacing was reduced progressively by pa 1 / v2 so that in each reduction the estimated leading finite 
lattice spacing artifacts were a factor of two smaller than in the previous set. The coarser lattices 
were added to support thermodynamics studies and to provide further leverage for continuum ex- 
trapolations. The medium coarse ensemble was added after coarse and fine and has a better tuned 
strange quark mass based on analysis of the other ensembles. 

For comparison, at a ~ 0.12 fm, a pa 0.09 fm and a ~ 0.06 fm, quenched ensembles with the 



10 In principle, the negative mass region can be simulated by adding a term to the action. Because of the sign 
problem, this wou ld be extremely challenging in four dimensions. However, it has been shown to work well in the 



Schwinger model dDiirr and Hoelblina 120061) 
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same gauge action were also generated. For each of these lattice spacings, the gauge coupling |3 = 
10/ g 2 was adjusted as the light quark mass was changed to keep the lattice spacing approximately 
fixed. However, the lattice spacing could only be determined accurately after the large ensembles 
were generated, so it is necessary to take into account the small differences in lattice spacing among 
the ensembles in the same set. In Sec. IIV.BI we then describe measurement of the lattice spacing 
on each ensemble, and a parameterized fit to smooth out statistical fluctuations. 

The strange quark mass in lattice units, am s , was estimated before simulations began, and was 
held fixed as the light quark mass and gauge coupling were varied. Later analysis, described in 
Sec. |VH determined the correct strange quark mass much more accurately, and in fact the initial 
estimates turned out to be wrong by as much as 25%. 

In the a ~ 0.12 fm set, several ensembles have a large dynamical quark mass — as large as 
eleven times the physical strange quark mass. This was done to investigate the physics of contin- 
uously "turning on" the quarks by lowering their masses from infinity. There are also a number 
of ensembles with a lighter-than-physical strange quark mass. These were generated to explicitly 
study dependence on the sea strange quark mass, and, since the lighter strange quark implies less 
sensitivity to higher orders in SU(3) chiral perturbation theory, enable improved determinations of 
the parameters in the chiral expansion, particularly of the low-energy constants (see SeclVTI). 

The fields satisfy periodic boundary conditions in the space directions, while the boundary 
condition in the Euclidean time direction is periodic for the gauge fields and antiperiodic for the 
quark fields. 

Table U shows the parameters of the asqtad ensembles (a few short "tuning" ensembles are not 
included). Here am/ is the dynamical light quark mass in lattice units and am s is the strange quark 
mass. Figure [8] plots the quark masses and lattice spacings of these ensembles. 



A. Algorithms and algorithm tests 



1987 ) de 



The earlier ensembles were generated using the "R" algorithm (jGottlieb et al. . 
scribed in Sec. III.Cl The molecular dynamics step size was generally set at about two thirds 
of the light quark mass in lattice units. More recent lattice generation has used rational func- 
tion approximations for the fractional powers described in Sec. III.Cl In those simulations, 
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P=10/g 2 


ami 


am s 


{L/af x (T/a) 


Lats. 


Uq 


r\/a 


m n L 


a « 0.18 fm 


6.503 


0.0492 


0.0820 


16 3 x48 


250 


0.85636 


1.778(8) 


9.07 


6.485 


0.0328 


0.0820 


16 3 x48 


334 


0.85585 


1.785(7) 


7.47 


6.467 


0.0164 


0.0820 


16 3 x48 


416 


0.85492 


1.801(8) 


5.36 


6.458 


0.0082 


0.0820 


16 3 x48 


484 


0.85489 


1.813(8) 


3.84 


a « 0.15 fin 


6.628 


0.0484 


0.0484 


16 3 x 48 


621 


0.8623 


2.124(6) 


8.48 


6.600 


0.0290 


0.0484 


16 3 x 48 


596 


0.8614 


2.129(5) 


6.63 


6.586 


0.0194 


0.0484 


16 3 x 48 


640 


0.8609 


2.138(4) 


5.46 


6.572 


0.0097 


0.0484 


16 3 x48 


631 


0.8604 


2.152(5) 


3.93 


6.566 


0.00484 


0.0484 


20 3 x 48 


603 


0.8602 


2.162(5) 


3.50 


a «0.12fm 


8.000 


OO 


CO 


20 3 x 64 


408 


0.8879 


2.663(6)* 


na 


7.350 


0.4000 


0.4000 


20 3 x 64 


332 


0.8822 


2.661(7)* 


29.4 


7.150 


0.2000 


0.2000 


20 3 x 64 


341 


0.8787 


2.703(7)* 


19.6 


6.960 


0.1000 


0.1000 


20 3 x 64 


340 


0.8739 


2.687(0)* 


13.7 


6.850 


0.0500 


0.0500 


20 3 x 64 


425 


0.8707 


2.686(8) 


9.70 


6.830 


0.0400 


0.0500 


20 3 x 64 


351 


0.8702 


2.664(5) 


8.70 


6.810 


0.0300 


0.0500 


20 3 x 64 


564 


0.8696 


2.650(4) 


7.56 


6.790 


0.0200 


0.0500 


20 3 x 64 


1758 


0.8688 


2.644(3) 


6.22 


6.760 


0.0100 


0.0500 


20 3 x 64 


2023 


0.8677 


2.618(3) 


4.48 


6.760 


0.0100 


0.0500 


28 3 x 64 


275 


0.8677 


2.618(3) 


6.27 


6.760 


0.0070 


0.0500 


20 3 x 64 


1852 


0.8678 


2.635(3) 


3.78 


6.760 


0.0050 


0.0500 


24 3 x 64 


1802 


0.8678 


2.647(3) 


3.84 


6.790 


0.0300 


0.0300 


20 3 x 64 


367 


0.8689 


2.650(7) 


7.56 


6.750 


0.0100 


0.0300 


20 3 x 64 


357 


0.8675 


2.658(3) 


4.48 


6.715 


0.0050 


0.0050 


32 3 x 64 


701 


0.8671 


2.697(5) 


5.15 



TABLE I Table of asqtad ensembles, uq is the input tadpole factor, Eq. Q, rather than the value determined 
from the ensemble average of the plaquette. Lattice spacings are from the smoothed fit described in the text, 
except where indicated by a "*". For these ensernl^es, r\/a is from this ensemble alone, rather than the 
smoothed fit. To convert to physical units, use r\ 0.31 fm. A f indicates that the run is in progress. This 
list of ensembles and counts of archived lattices are as of December 2008. 
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ami 


am s 


(L/a) 3 x (T/a) 


Lats. 


Uq 


n/a 


m n L 


a w 0.09 fm 


8.400 


OO 


OO 


28 3 x 96 


396 


0.89741 


3 730(7)* 


na 


7.180 


0.0310 


0.0310 


28 3 x 96 


500 


8808 


3 822(10) 


8.96 


7.110 


0.0124 


0.0310 


28 3 x 96 


1996 


0.8788 


3.712(4) 


5.78 


7.100 


0.0093 


0.0310 


28 3 x 96 


1138 


0.8785 


3 705(3) 


5.04 


7.090 


0062 


0.0310 


28 3 x 96 


1946 


0.8782 


3.699(3) 


4.14 


7 08 s ) 


00465 


0310 


32 3 x 96 


540 1 " 


8781 


3 697(3) 


4.11 


7.080 


0.0031 


0.0310 


40 3 x 96 


1012 


0.8779 


3.695(4) 


4.21 


7.075 


0.00155 


0.0310 


64 3 x 96 


530 f 


0.877805 


3.691(4) 


4.80 


7.100 


0.0062 


0.0186 


28 3 x 96 


985 


0.8785 


3.801(4) 


4.09 


7.060 


0.0031 


0.0186 


40 3 x 96 


642 


0.8774 


3.697(4) 


4.22 


7.045 


0.0031 


0.0031 


40 3 x 96 


440 t 


0.8770 


3.742(8) 


4.20 


a » 0.06 fm 


7.480 


0.0072 


0.0180 


48 3 x 144 


625 


0.8881 


5.283(8) 


6.33 


7.475 


0.0054 


0.0180 


48 3 x 144 


617 


0.88800 


5.289(7) 


5.48 


7.470 


0.0036 


0.0180 


48 3 x 144 


771 


0.88788 


5.296(7) 


4.49 


7.465 


0.0025 


0.0180 


56 3 x 144 


800 


0.88776 


5.292(7) 


4.39 


7.460 


0.0018 


0.0180 


64 3 x 144 


826 


0.88764 


5.281(8) 


4.27 


7.460 


0.0036 


0.0108 


64 3 x 144 


483 


0.88765 


5.321(9) 


5.96 


a « 0.045 fm 


7.810 


0.0028 


0.0140 


64 3 x 192 


861 


0.89511 


7.115(20) 


4.56 



TABLE II Table Ucontinued. 



we have used the Omelyan second order integration a 



Sexton and Weingartenl 



1992 



the fermion and gauge forces (|Sexton and Weingarten , 



akaishi and de Forcrand 



gorithm (lOmelyan et all 



2002a b 



2003 



2006). We used different step sizes for 



19921) . with the step size for the fermion 



force three times that of the g auge force. We used four sets of pseudofermion fields and cor 



responding rational functions (|Hasenbusch . 



2001 



Hasenbusch and Jansenl . 



20031) . The first set 



implements the ratio of the roots of the determinants for the physical light and strange quarks to 
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FIG. 8 Lattice spacings and quark masses used. The octagons indicate ensembles with the strange quark 
near its physical value, while the crosses indicate those with an unphysically light strange quark. The burst 
at lower left shows the physical light quark mass. Here the quark masses are in units of MeV, but using the 
asqtad action lattice regularization. 

the determinant for three heavy "regulator" quarks with mass am r = 0.2. That is, it corresponds 
to the weight det(M(m/)) 1 / 2 det(M(m i y)) 1 / 4 det(M(m r )) _3 / 4 . The remaining three pseudofermion 
fields each implement the force from one flavor of the regulator quark, or the fourth root of the cor- 
responding determinant. These choices are known to be reasonably good, but could be optimized 
further. 

For all but the largest lattices generated with rational function methods, we included the 
Metropolis accept/reject decision to eliminate step size errors, or the RHMC algorithm. Because 
the integration error is extensive, use of the RHMC algorithm for the largest lattices would have 
forced us to use very small step sizes and double precision in many parts of the integration. For 
these lattices it was much more efficient to run at a small enough step size that the integration error 
was less than other expected errors in the calculation (the RHMD algorithm). 
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FIG. 9 The plaquette as a function of integration step size squared for 20 3 x 64 lattices with P = 6.76 and 
am q = 0.01/0.05. The point at £ 2 = is from the RHMC algorithm, and the point indicated by R is the 
v alue used in the R alg orithm production runs. The remaining two points are from short test runs described 



in 



Aubin et al. 



(2004a). 



Errors from the integration step size in the R algorithm were o rigin ally estimated from s hort 



runs with different step sizes, as described in 



Bernard et al. 



( 2001 ) and 



Aubin et al. 



(l2004al) . In 



several cases, ensembles originally generated with the R algorithm were later extended with the 
RHMC algorithm. This allows an ex post facto test of the step size errors in the R algorithm, with 
much higher statistics than possible for a tuning run. Figure [9] shows the average plaquette for one 
a ~ 0.12 fm run as a function of step size squared, combining the early tuning runs with the R and 
RHMC algorithm production runs. Table HITl compares the expectation values of the plaquette and 
the light quark condensate and, in some cases, the lattice spacing and pion mass, for the ensembles 
where both algorithms were used. The differences are small and in most cases are comparable to 
the statistical errors. 

In one case, a « 0. 12 fm and am q = 0.01/0.05, an ensemble with a larger spatial size (28 3 ), was 
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0.020 
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052553(611 052306(281 00025 1(671 


6.76 


0.010 


0.050 


0.00667 


1.700917(21) 1.700879(18) 


-0.000038(28) 


0.036875(43) 0.037174(36) 0.000300(56) 


6.76 


0.007 


0.050 


0.00500 


1.701183(22) 1.701177(18) 


-0.000006(29) 


0.031388(54) 0.031306(38) -0.000082(66) 


6.76 


0.005 


0.050 


0.00300 


1.701181(17) 1.701211(11) 0.000030(20) 


0.027551(50) 0.027597(25) 0.00045(56) 


7.11 


0.0124 


0.031 


0.00800 


1.789213(19) 1.789075(7) 


-0.000138(20) 


0.024584(22) 0.024620(10) 0.000036(24) 


7.09 


0.0062 


0.031 


0.00400 


1.784552(9) 


1.784541(6) 


-0.000011(11) 


0.015622(17) 0.015608(14) -0.00015(22) 
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0.00200 


1.782300(8) 


1.782254(11) 


-0.000046(11) 


0.010664(18) 0.010860(19) 0.000196(26) 
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^ (RHMC) 


difference 


am n (R) am % (RHMC) difference 


7.11 


0.0124 


0.031 


0.00800 


3.708(13) 


3.684(17) 


-0.024(21) 


0.20640(20) 0.20648(20) 0.00008(28) 


7.09 


0.0062 


0.031 


0.00400 


3.684(12) 


3.681(8) 


-0.003(14) 


0.14797(20) 0.14767(13) -0.00030(24) 


7.08 


0.0031 


0.031 


0.00200 


3.702(8) 


3.682(7) 


-0.020(11) 


0.10528(9) 0.10545(9) 0.00017(13) 



TABLE III Comparison of plaquette and light quark condensate for ensembles run partly with the R algo- 
rithm and partly with the RHMC algorithm. For the a « 0.09 fm ensembles, we also show r\/a and the pion 
mass. 



generated to check for effects of the spatial size. In general, these effects were found to be small as 
expected, although the effects on f % and fx differ significantly from one-loop chiral perturbation 
theory estimates, as discussed in Sec.lVTl 



B. The static potential and determining the lattice spacing 

Since results of lattice QCD simulations are initially in units of the lattice spacing, knowing 
the lattice spacing is crucial to calculating any dimensionful quantity. However, since ratios of 
dimensionful quantities (mass ratios) calculated on the lattice will only have their physical values 
at the physical quark masses and in the continuum limit, there is arbitrariness in the determination 
of the lattice spacing except in the physical limit. Some dimensionful quantity must be taken to be 
equal to its physical value or to some a priori model. 

Fol lowing the pra ctice of most current lattice simulation programs, we use a Sommer 



scale (|Sommerj . 



19941) as the quantity kept fixed, and determine this scale from some well con- 
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FIG. 10 The static quark potential for the ensemble with a « 0.09 fm and m; 0.2m s . This was obtained 
from time range five to six. The inset magnifies the short distance part, showing a lattice artifact which is 
discussed in the text. 

trolled measurement. 

A Sommer scale is defined as the length where the force between a static (infinitely heavy) 
quark and antiquark satisfies r 2 F(r) = C, where C is a constant. Intuitively, this is a length where 
this static potential changes behavior from the short distance Coulomb form to the long distance 
linear form. In particular, the most common choice is rpj, defined by C = 1.65. We have chosen to 
use r\, defined by C = 1. This choice was made based o n early simula t ions at a ~ 0.12 fm where 



it was found that r\ had smaller statistical errors than r$ (|Bernard et al. 



2000b) 



The calculation of the static potential on the earlier ensembles is described in 



Bernard et al. 



(I2000b|) . We begin by fixing the lattice to Coulomb gauge. In this gauge, we can evaluate the 



potential from correlators of (nonperiodic) Wilson lines, where the line at (x,t) with length T is 
Wj{x,t) = YiJ^Q 1 U4(x,t + i). The Coulomb gauge fixing, which makes the spatial links as smooth 
as possible, is an implicit way of averaging over all spatial paths closing the loop at the top and 
bottom. Because we do not explicitly construct the spatial parts, it is easy to average over all lattice 
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FIG. 1 1 The static quark potential and first excited state potential for the ensemble with a « 0.06 fm and 
mi ~ 0.lm s . This was obtained from time range three to twenty, using the APE smeared time links discussed 
in the text. 

points (x,t) and to get the potential at all spatial separations R. 

The first step in determining r\ is to extract V(R) from the expectation value of the correlators 
of Wilson lines. We expect 

L{RJ) = (w}(x,t)W T (x + R,t)^ = Ae- V{R)T +A'e- V '^ )T + ... , (113) 

where V', etc. are potentials for excited states. For a > 0.09 fm, the excited states are negligible 
for fairly small T, and we simply take V(R) = \og(L(R,T) /L(R,T + 1)). Specifically, we use 
T = 3 for a « 0. 15 fm, T = 4 for a w 0. 12 fm and T = 5 for a w 0.09 fm. Figure Ql shows the 
resulting potential for the run at a ~ 0.09 fm and m/ = 0.2m s . The inset in this figure shows the 
short distance part of the potential. In this inset, there is a visible lattice artifact where the point 
at R = 2, or separation (2,0,0) is slightly below a smooth curve through the nearby points with 
off-axis distances R. However, at R = 3 the lattice artifacts are quite small. In fact, what appears to 
be a single point at R = 3 is actually two points, one for R = (3, 0, 0) and another for R = (2,2,1). 
The small objects in the center of the plot symbols are the statistical error bars on V(R). 
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FIG. 12 The static quark potential in units of r\ for five different lattice spacings. In all cases, these are for 
light quark mass of two tenths the simulation strange quark mass. For each lattice spacing, a constant has 
been subtracted to set V(r\) = 0. The ruler near the bottom of the plot shows distance in units of fm, using 
r\ = 0.318 fm. The multiple rulers in the upper half of the plot show distance in units of the lattice spacings 
for the different ensembles. 

For a ~ 0.06 fm, the above procedure for finding V(R) gave large statistical errors. This is 
primarily because a large constant term in the potential causes a rapid falloff of L(R,T) with T. 
This constant can be considered a self energy of the static quark, diverging as I /a. To fix this, th e 



timelike links were smeared by adding a multiple of the three link "staples" (lAlbanese et al. 



,11282), 



namely "fat3 links" defined in Eq. (|69l ) with G) = 0.1. The Wilson line correlators L{R, T) were 
computed from the smeared time direction links as described above. As expected, this reduces the 
constant term in V(R), and comparison with the potential from unsmeared links suggests that any 
systematic effects on r\/a are less than 0.005 at a ~ 0.06 fm, smaller than the statistical errors. 
With the smeared time links, the correlators L(R 1 T) are statistically significant out to T as large 
as twenty (for small R). It is then advantageous to do a two state fit to L(R, T). For the a ~ 0.06 
fm ensembles, we generally fit these two states over a time range 3 < T < 20. An example of the 
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potential from this procedure is shown in FigQTJ The first excited state potential is also shown, but 
we caution the reader that in addition to having large statistical errors this excited state potential 
has not been carefully checked for stability under varying fit ranges, or under addition of a third 
state to the fit. 

Once V(R) is determined, we find r\ by fitting V(R) to a range of R approximately centered at 
r\. We use a fit form 

V(*)=C+| + aR + xQ -£) (114) 

Here C is part of the quarks' self energy, o is the string tension and B is ^(X s for a potential 
definition of a s . The last term, t — ^, is the difference between the lattice Coulomb potential, 
Tt\ lot = / (f^^oo (p) e ' pR w i m £*oo (?) me f ree lattice gluon propagator calculated with the 
Symanzik improved gauge action, and the continuum C oulomb potential l/ R. Use of this correc- 



tion term was introduced by the UKQCD collaboration (jBooth et at 



19920 . This correction was 



used for R < 3. The scale r\ (or tq) was then found from solving r 2 F(r) = C with X set to zero, 
r\ = J ii^. Since we often want lattice spacing estimates from only a few lattices, and there are 
a large number of distances to be fit, these fits were generally done without including correlations 
among the different R. Errors on r\ are estimated by the jackknife method, where the size of the 
blocks eliminated ranges from 30 to 100 simulation time units. Spot checks comparing fits includ- 
ing the correlations confirmed that the jackknife errors are consistent with derivative errors in the 
correlated fits, and that the fit function does fit the data well over the chosen range. 

For the a ~ 0. 1 8 fm ensembles, we used the spatial range from 1 .4 or 1.5 to 6.0; for the a ~ 0. 15 
fm ensembles, y/2 <R<5; for the a ~ 0. 12 fm ensembles, y2 <R<6; and for the a « 0.09 fm 
ensembles, 2 < R < 7. For the a ~ 0.06 fm ensembles, where the two state fits with smeared links 
were used, the spatial range was 4 < R < 7, and for the a ps 0.045 fm run, it was 5 < R < 10. 

The static quark potentials for different lattice spacings can be overlaid after rescaling to check 
for lattice effects and to plot the potential over a large range. Figure [[2] shows such a plot in 
units of r\ for five different lattice spacings, using the ensembles with m\ = 0.2m s at each lattice 



spacing. In lBernard et al.\ (|2000br) . it was found that including the dynamical quarks modifies the 
static potential in the expected way. This can be seen by plotting dimensionless quantities such as 
ro / r\ or r\ y/o. When this is done in a region where the potential is approximated by Eq. (II 14|) and 



r\ is found from r\ = y -^p, this amounts to plotting the coefficient of l/R in the fit. 

Once r\ is estimated for each ensemble, the estimate can be improved by fitting all values of 
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r\/a to a smooth function of the gauge coupling and quark masses. We have used two different 
forms for this smoothing. In the first form, we fit logfn /a) to a polynomial in [3 and lanii + am s . 



The second form is a function based on work of 



AUton 



(1996): 



a_ = Cof + C 2 g 2 f 3 +C 4 g 4 f 



where 



f = (b g 



2N(-*!/(2fc§)) 



exp(-l/(2%> 2 )) , b = (11 -2n//3)/(4jc) 



(115) 



bi = (102-38/z//3)/(47c) 4 
Q = C o + Couami/f + Coi s am s /f + Co 2 (am tot ) 2 



am tot = lami / f + am s / f 
C 2 =C 2 o + C 2 iam tot . 



(116) 



Here Coo, Ql/, Q)i.v, C02, C20, C21, C4, and D 2 are parameters. The second form is a slightly better 
fit, and we have used it for the r\/a values in TablelH Errors on the smoothed r\ ja are estimated by 
a bootstrap for which artificial data sets were generated. In these data sets the value ofr\/a for each 
ensemble was chosen from a Gaussian distribution centered at the value for the ensemble given by 
the fit, and the standard deviation was given by the statistical error in r\ja for the ensemble. 

To find r\ in physical units, we use a quantity that is both well known experimentally and ac- 
curately determined in a lattice calculation. One such quantity, and the one used in most of our 
work, is the splitting between two energy levels of the bb mesons. These splittings ha ve been cal 



culated on several of the a s qtad e nsembles by the HPQCD/UKQCD collaboration (IGray et all 



2003, 



2005; 



Wingate et all 120041) . From fitting the 2S-1S splittings on the a w 0.12 fm en- 



sembles with quark masses ami/am s = 0.01/0.05, 0.02/0.05, 0.03/0.05 and 0.05/0.05, and the 
a w 0.09 fm ensembles with light masses ami/am s = 0.0062/0.031 and 0.0124/0.031, to the form 
r](a,ami,am s ) = r\ ys + c\a 2 + c 2 ani[/am s , we find r^ hys = 0.318 fm with an error of 0.007 fm. 



( Gray et al. 



( 2005 ) used a different fitting procedure to estimate r\ hys = 0.321(5) fm.) 
More recently, analysis of the light pseudoscalar meson masses and decay constants gave 
an accurate value of f % . The fitting procedure to arrive at this is complicated — see Sec. |Vll 
Requiring that f n in the continuum and chiral limits match its experimental value gives r\ = 
0.3108(15)(±^) fm, where the errors are statistical and systematic, respectively. 

To summarize, we set the scale for each ensemble by a = {a/r\) x rf* ys , where (a/r{) is the 
output of the smoothing function, Eq. (II 151) . at the ensemble values of ami, am s , and g 2 , and r^ hys is 
the physical value of r\, obtained either from bb mesons splittings or f n . The scheme is useful for 
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generic chiral extrapolations, and tends to result in fairly small dependence of physical quantities 
on the sea-quark masses. However, chiral perturbation theory assumes a mass -independent scale 
setting scheme, because all dependence on quark masses is supposed to be explicit. So detailed 
fits to chiral perturbation theory forms require a mass -independent scale procedure, especially if 
one hopes to extract low energy constants that govern mass dependence. Once the r\ smoothing 
form is known, though, it is easy to modify the procedure to make it mass independent: instead of 
using the ensembles' values of am/ and am s in the smoothing function, Eq. ( II 151 ), use the physical 
values. This mass-independent scheme is used for the analysis of light pseudoscalars described in 
Sec. El 



C. Tuning the strange quark mass 

In most of these ensembles, the original intent was to fix the strange quark mass at its correct 
value, and to set the light quark mass to a fixed fraction of the strange quark mass. The correct 
strange quark mass, however, is actually not known until the lattices are analyzed. In practice, the 
best that can be done is to estimate the correct strange quark mass from short tuning runs or by 
scaling arguments from results of earlier runs. As described in Sec. |VlJ the physical strange and 
up/down quark masses are determined by demanding that the light pseudoscalar meson masses 
take their physical values. For the strange mass, we find am s = 0.0439(18) at a ~ 0. 15 fm, am s = 
0.0350(7) at a « 0.12 fm, am s = 0.0261(5) at a « 0.09 fm and am s = 0.0186(4) at a « 0.06 fm. 
For the up/down mass, we find am/ = 0.00158(7) at a ~ 0.15 fm, am/ = 0.00126(2) at a ~ 0.12 
fm, ami = 0.000955(8) at a « 0.09 fm and am, = 0.000684(8) at a « 0.06 fm. The errors include 
statistical and systematic effects, but they are dominated by the systematic effects. 



D. The topological susceptibility 

The topological structure of the QCD vacuum is an important characteristic of the theory. A 
stringent test for lattice simulations consists in correctly capturing the dependence of the topo- 
logical susceptibility on the number of quarks and their masses, since this susceptibility reveals 
the effect of the quarks on the nonp erturbative vacuum structure. Chiral perturbation theory pre- 



dicts %topo(nf,mi) in the chiral limit (ILeutwyler and Smilga . 



19921) . Lattice calculations, however, 



have struggled to reproduce this dependence satisfactorily because the topological charge is not 
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uniquely defined and the fermion action typically breaks chiral symmetry. The asqtad action com- 
bined with rSZPT gives us good control over the taste and chiral symmetry breaking effects; thus 
we expect that a careful treatment of the topological charge will lead to an accura t e computation 



of the topo 



(2004 and 



ogical suscep t ibility. This has been explored in 



Bernarc 



As explained in 



et al. 



Bernard et al. 



(20071) 



(I2003dh . 



Billeter et al. 



Billeter et al. 



(|2004|) . the chiral anomaly cou- 



Aubin and Bernard! (12003 al) and 
pies to the taste-singlet meson, not the Goldstone pion, which is the usual focus of hadron spec- 
troscopy calculations. (Of course, in the continuum limit these mesons are degenerate.) To leading 
order in rSZPT, the topological susceptibility depends on this mass as 

%top ° " 1+^/(^+3^/(2*® ' ° ^ 

where m n j is the taste-singlet pion mass, and mo comes from the term representing the coupling 
of the anomaly to the r\' in the chiral Lagrangian, Eq. (l93l) . The strange flavor-singlet, taste-singlet 
meson mass is denoted m ss j. 
E quation (Qj 



tion (|Venezianol . 



71) interpolates smoot hly between the infinite sea-quark-mass (quenched) predic- 



1979 



Witten . 



19791) . % = f^m^/Yl, which we can use to set mo, and the chiral 
limit, m/ — > 0, which is dominated by the pion, X = f^rn^/S. Hence, to this order, we simply 
replace the Goldstone pion mass with the mass of the taste-singlet (non-Goldstone) pion in the 
Leutwyler-Smilga formula. Note that this means that, at nonzero lattice spacing, the topological 
susceptibility does not vanish as m/ — > 0, a reminder that the continuum limit must be taken before 
the mi — ► extrapolation. 

In order to compute the topological charge den sity q(x) on our lattice ensembles, we use three 



iterat ions of the Boulder HYP smoothing me t hod (IDeGrand et al. . 



de Forcrand et al 



20011) . wh ich we have found ([Bernard et al. . 
method of 



1997 



Hasenfratz and Knec htli. 



2003alidl) compares well with the improved cooling 



(|1997|) . We define the topological susceptibility from the correlator 



of q(x) via 



Xto P o = (£ 2 )/V = J d 4 r(q(r)q(0)) . 



(118) 



On our lattices, the short-distance part of the density correlator has a strong signal, but the cor- 
relator at large separation is noisy. To reduce the resulting variance, we define a cutoff distance 
r c . In the integral above, for r < r c where the signal is strong, we use the measured values of the 
correlator (q(r)q(0)). For r > r c we integrate a function obtained by fitting the measured correlator 
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FIG. 13 Points used to compute (q(r)q(O)). Measured points (open symbols') are u sed for r < r c ~ 9a. For 



r > r c the fit function (solid curve) is used in Eq. (1118b - From 
to a Euclidean scalar propagator 



Bernard et al. 



(I2007fh . 



(q(r)q(0) ) ~ A^Ki (m^r) /r +A r] ,K l (m^r) jr 



(119) 



where we use priors for the masses of the T| and r\', and K\ is a Bessel function. This significantly 
reduces the variance in Q 2 . An example of the measured values of q(r), the fit function, and the 
fitting range are shown in Fig. IT3l 

Figure Qj] shows this definition of Xtopo computed on our coarse (a ~ 0.12 fm), fine (a ~ 0.09 
fm), and superfine (a ~ 0.06 fm) lattices. The continuum limit is taken first by fitting the suscepti- 
bility data to 



1 



%topo o 



[m 



A +Aia 2 + (A 2 +A 3 a 2 +A 4 a 4 )/m^ / . 



(120) 



The solid black line in Fig. [14] shows the a — ► form of this function. Some representative 
points along this line are shown with error bars reflecting the errors of the continuum extrapolation. 
Finally, the chiral perturbation theory prediction of Eq. (II 171 ), shown as a dotted line, is based on 
the value for mo set by the quenched data. 

With the addition of the new a ~ 0.06 fm data, we see that the topological susceptibility is 
behaving as expected in the m 2 j^0 limit of rooted staggered chiral perturbation theory. 

These results lend further credibility to the use of the fourth root procedure to simulate single 
flavors, since aberrant results from this procedure would be expected to arise first in anomalous 
behavior of topological quantities and correlations, as these are rather sensitive to the number of 
flavors. 
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FIG. 14 Topological susceptibil ity data, and it s contin uum extrapolation, compared with the prediction of 



Eq. <1 1 17b . Update of figure from 



Bernard et al. 



<|2007j). 



V. SPECTROSCOPY OF LIGHT HADRONS 



Computing the masses of the light hadrons is a classic problem for lattice QCD, since the 
masses and structures of these particles are highly nonperturbative. By this point, hadron mass 
computations, including the effects of light and strange dynamical quar ks, have been done for 



several different lattice a ctions, including staggered quarks, Wilson quarks (IDiirr et al. 



Ukita et all 



2007 



2009|) and domain- wall quarks (|Allton et al. . 



2008 



Ukita et al. 



2008, 



2009; 



2007). It has 



long been apparent from these and other studies that lattice QCD reproduces the experimental 
masses within the accuracy of the computations. For most of the light hadrons, however, this 
accuracy is not as good as for many of the other quantities discussed in this review. The reasons for 
this are that these masses have a complicated dependence on the light quark mass, making the chiral 
extrapolation (to the physical light quark mass) difficult, and that all but a few of these hadrons 
decay strongly. Most of the lattice simulations are at heavy enough quark masses or small enough 
volumes that these decays cannot happen, so the chiral extrapolation crosses thresholds. With 
staggered quarks there is the additional technical complication that for all but the pseudoscalar 
particles with equal mass quarks the lattice correlators contain states with both parities, with one 
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of the parities contributing a correlator that oscillates in time. 

Masses of the lowest-lying light-quark hadrons have been computed on almost al l of the MILC 



asqtad ensembles. Hadron masses from the a ~ 0. 12 fm ensem bles were reported inlBernard et al. 



(1200 lb . masses from the a ps 0.09 fm ensembles w ere added in 



Aubin et al. 



Bernard et al 



(I2004a|) . and nucleon 



(|2007cb . Simple extrapolations 



and Q. masses from the a ~ 0.06 fm ensembles in 
of these masses to the continuum limit and physical quark mass, including results from several of 
the a ps 0.06 fm ensembles, are compare d to experiment in Fig. 



charm and bottom meson mass splittings (IGrav et al. , 



with experimental values (|Amsler et al. 



2003. 



2005 



5J In addition, this fig ure shows 



Wingate et all 120041) compared 



2008). 



A. Hadron mass computations 



The theory behi n d hadr o n mass computation s wit h staggered quarks was de veloped in 



Klub erg-Stern et al. 



Kilcup and Sharpe ( 



elude 



(Il983al) . 



Marinari et al. 



Goltermanl (I1986bl) and 



Golterman and Smitl (11985b (see also 



1987b). E arly implementation s ., in which te c 



1981ab . 



Bowler et al. 



(Il987b . 



Gupta et al. 



inical aspe cts were addressed, in 



(11991b . and 



Fukugita et al. 



(1 1993b . 



The calculation begins with a Euclidean-time correlation function for any operator that can 
produce the desired state from the vacuum. For instance, if an operator O can annihilate a particle 
p and the adjoint can create p, then we study the zero-momentum correlation function, or 
"correlator" C t given by 

C oto (?)=£(o(x,0o t (0,0)). (121) 

X 

By putting in a complete set of states between the two operators, we find 



c o t o (0 



£(0| O \n) (n\ O f |0) exp(-M„t) 



(122) 



If the particle p is the lowest-energy state n, then for large Euclidean time, the dominant contribu- 
tion will be |(0|o|p)| 2 exp(— M p t). Generally, there will be additional contributions from higher 
mass states, and with staggered quarks there are usually contributions from opposite parity states 
of the form (— l) f exp(— M't). In addition, because of the antiperiodic boundary conditions in time 
for the quarks, there will be additional terms of the form exp(— M n (T — t)), where T is the time 
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FIG. 15 The "big picture" — comparison of masses calculated on the asqtad ensembles with experimental 
values. For the light quark hadrons we plot the hadron mass, and for the cc and bb masses the difference 
from the ground state (IS) mass. The continuum and chiral extrapolations of the pion and kaon masses are 
described in Sec. |VIJ and most other meson masses were extrapolated to the continuum and physical light 
quark masses using simple polynomials. Masses of hadrons containing strange quarks were adjusted for the 
difference in the strange quark m ass used in generating the ensembles from the correct value. The nucleon 



mass extrapolation, described in 



Bernard et al. 



The charmonium ma ss spl itting is fro m 



Wingate et al. 



(2004)) and 



Gray et al. 



Follana et al. 



(2007 d), use d a one-loop chiral pertur bation theory form . 



d2008l) . and the bb split tings from 



Gray et al. 



Amsler et al. 



(2003). 



(120081) . TheT 



(2005). Experimental values are from 
2S-1S splitting and the % and K masses are shown with a different symbol since these quantities were used 
to fix r\ in physical units and the light and strange quark mass es. Earlier versions of the plot appeared in 



Aubin et al. 



(I2004al) and the PDG "Review of Particle Physics" (lAmsler et al. 



2008). 
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FIG. 16 Pion and nucleon correlators plotted vs. the distance from the source. These correlators are from 
the P = 6.76, ami/am s = 0.007/0.05 ensemble. The small symbols in the center of the octagons in the pion 
correlator are error bars. Note the increasing fractional errors with distance in the nucleon correlator, and 
the constant fractional errors in the pion correlator. 

extent of the lattice. Thus, with staggered quarks a meson correlator generically has the form 

C oto (0 = Ao{e- M ^ + e- M ^ T - t ^+A^e- M ^ + e- M ^ T - r ^+... 

+ (-l)%( e - M o t + e- M 'o( T - t A+... (123) 

Here the primed masses and amplitudes with the factor of (— l) r correspond to particles with 
parity opposite that of the unprimed. For baryons the form is similar, except that the backwards 
propagating terms {e~ MtyT ~'^) have an additional factor of (— l) f+1 . Here the overall minus sign 
in the backwards propagating part is due to the antiperiodic boundary conditions for the quarks in 
the Euclidean time direction. Figure [T6] shows correlators for the pion and nucleon in a sample 
asqtad ensemble. Statistical errors on the pion correlator are the tiny symbols in the center of the 
octagons. The effect of periodic (for a meson correlator) boundary conditions in time is clearly 
visible. For short times, there are contributions from heavier particles. 

For hadrons other than glueballs, evaluating this correlator requires computing M~y where M 
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is the matrix defining the quark action. This can be done by making a "source" vector b which 
is nonzero only at lattice point y, or in some small region, and solving the sparse matrix equa- 
tion Ma = b, usually using the conjugate gradient algorithm. (Here a and b are vectors with one 
component for each color at each lattice site in the system - i.e., 3V complex components. With 
Wilson-type quarks there would also be four spin components per lattice site.) 

The simplest possibility for is an operator built from quarks and antiquarks located in the 
same 2 4 hypercube, often even on the same lattice site. This is usually called a point source. 
Because the point operator Op tends to have a large overlap with excited states, it is usually advan- 
tageous to take a "smeared" source operator o \ where the quarks in the hadron may be created at 
different lattice sites. One common approach is to choose a smeared operator that creates quarks 
and antiquarks with a distribution similar to that of the expected quark model wave function of the 
desired hadron. A cruder and simpler approach used in most of the MILC light hadron mass calcu- 
lations is to take a "Coulomb wall" source, where the lattice is first gauge transformed to the lattice 
Coulomb gauge, making the spatial links as smooth as possible. Then a source is constructed 
which covers an entire time slice, for example, with a 1 in some corner of each 2 3 cube in the 
time slice. This works because with Coulomb gauge fixing contributions from source components 
within a typical hadronic correlation length interfere coherently, while contributions average to 
zero if the quarks created by are widely separated (although they do contribute to the statistical 
noise). In other words, (m^\ ^ t MZ^ 2 f y is significant only when \x\ — X2\ is less than a typical 
hadronic size. For example, a Coulomb wall operator appropriate for a Goldstone pion is 

M0 = Ex&0(-i)**x&0 • ( 124 ) 

x,y 

In a mass calculation, we want the state with zero spatial momentum, which is isolated by 
summing the sink position over all spatial points on a time slice. In many matrix element studies, 
we need hadron states with nonzero momenta, and they are isolated by summing over the spatial 
slice with the appropriate phase factors. 

Statistics are usually further enhanced by averaging correlators from wall sources, or other types 
of sources, from several time slices in the lattice. In general, each different source requires a new 
set of sparse matrix inversions. 

For most hadrons, statistical error is the limiting factor in the mass computations. At long 
Euclidean time t, a correlator with hadron H as its lowest mass constituent is proportional to 
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e Mlit . The variance of this correlator can itself be thought of as the correlator of the square of the 
operator 



H {x)0] i {x) 0l{y)0 H {y) 



(125) 



ladro ns it is understood that quark lines all run 



199(1) . The behavior of the variance at long dis- 



where in this correlator for flavor-no nsinglet 
from the operators at x to those at y (|Lepage . 
tances is dominated by the lowest mass set of particles created by Oh{x)O h (x). Thus for mesons 
Oh(x)oJ 1 (x) creates two quarks and two antiquarks which can propagate as two pseudoscalar 
mesons. Then the variance decreases approximately as e~ 2Mpst , where Mps is the mass of the 
pseudoscalar meson made from the quarks in 0^0 For baryons there are three quarks and three 
antiquarks, and the variance decreases approximately as e~ 3Mpst . This behavior can be seen in 
Fig. EH where the fractional error on the pion correlator does not increase with distance, while the 
fractional error on the nucleon correlator grows quickly. 

As discussed in Sec. III.B.3L hadrons with staggered quarks come with different "tastes," all 
of which are degenerate in the continuum limit. For pseudoscalar mesons, the mass differences 
between different tastes are large, but they are well understood as discussed in Sec. IIII.Al For 
the other hadrons, for which chiral symmetry is not the most important factor in determining the 
mass, taste symmetry violations are much smaller. In particular, we have computed masses for four 
different tastes of the p meson on many of our ensembles , and have failed to find any statistically 



significant taste splittings. (See also 



Ishizuka et al. 



(ll994lU 



B. Correlated fits 



There are several kinds of correlations in the numerical results of lattice gauge theory simu- 
lations. The Markov chain that produces the configurations produces correlated configurations. 
Thus, there are correlations in "simulation time." The correlations vary with the algorithm, and 
one can reduce them by increasing the simulation time gap between the configurations that are 
analyzed. Generation of configurations is computationally expensive, however, and the autocor- 
relation length is unknown until the run and some analysis is completed, so one usually saves 
configurations with some degree of correlation. A simple way to deal with these correlations is to 
block successive configurations together and then to estimate errors from the variance of blocks. 
However, if the number of blocks is not many times larger than the number of degrees of freedom, 
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the finiteness of the sample size must be considered when estimating goodness- of-fit or statistical 



errors on the parameters in a fit ([Michael . 



1994; 



Toussaint and Freeman 



2008). In cases where 



blocking is not practical, notably the pseudoscalar meson analysis in Sec. [VH we have estimated 
elements of the covariance matrix by using the measured autocorrelations in the data to rescale a 
covariance matrix based on unblocked data. 

Even if successive configurations are not correlated, different physical quantities are correlated 
with each other. For example, if the pion correlator is larger than average at a separation t from the 
source on a particular configuration, it is likely to be larger at t + 1. Thus, when extracting hadron 
masses, or other fit parameters, we must use the full correlation matrix in the fit model, not just the 
variance in each particular element fit. To be specific, let the values of the independent parameters 
be denoted xi and corresponding lattice "measured" values be yi. The fitting procedure requires 
varying the model parameters {A,} that define the model function yui^U {M) m order to minimize 
% 2 . For uncorrelated data, 

x 2 = %(yM{xum)-yi? J°h (126) 

where a, is the standard deviation of y ; . When the data is correlated, let Qj = Cov(y,, yj) and then 

X 2 = ~yi)Cr/ (y M (xj,{X}) - yj ) {111) 

hi 

(In practice C, 7 is almost always estimated from the same data as the y,, and in this case % 2 is more 
properly called T 2 .) Uncorrelated data reduces to Cu = 5 (J o 2 . If Cu has positive off-diagonal 
entries, then the data will look smoother than it would if uncorrelated. 

In Fig. \T7\ we show how the fitted pion and nucleon masses vary with the minimum distance 
from the source that is included in the fit. The octagons and squares are correlated fits, minimizing 
% 2 in Eq. (11271) . For the pion, the octagons correspond to a single-particle (two-parameter) fit, and 
the squares correspond to a two-particle (four-parameter) fit. For the nucleon, the octagons are fits 
including one particle of each parity. We need to decide which fit is best, and we do that based on 
the confidence levels of the fits, which are roughly indicated by the symbol size. Figure [FT] also 
contains fits ignoring correlations while minimizing the % 2 in Eq. (11261) . It can be seen that the error 
bars on these points are in general incorrect — they are neither a correct estimate of how much the 
parameters would likely vary if the calculation were repeated, nor of how much the parameters are 
likely to differ from the true value. We also see that the confidence levels are generally too large for 
the uncorrelated fits. In particular, based on its confidence level, one might accept the uncorrelated 
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FIG. 17 Result of fitting the correlators in Fig. [16] from a minimum distance to the center of the lattice 
(for the pion) or distance at which the correlator loses statistical significance (for the nucleon). For the 
pion correlator (left panel), octagons correspond to single-particle fits and squares to two-particle fits. The 
diamonds are from single-particle fits ignoring correlations among the data points. For the nucleon fits (right 
panel), all the fits use two particles, one of each parity. Octagons are correlated fits, and diamonds are fits 
ignoring the correlations. The sizes of the symbols are proportional to the confidence level of the fits, with 
the symbol size in the legends corresponding to 50% confidence. 

pion fit with minimum distance five. But in fact it can be seen that it differs significantly from 
the asymptotic value. The effects on the confidence level from ignoring correlations can be quite 
extreme. For example, in the single-particle pion fits with D min = 5, the correlated fit has % 2 = 180 
for 25 degrees of freedom, for a confidence of 10~ 24 , while the uncorrelated fit has % 2 = 14 for 25 
degrees of freedom, or an (erroneous) confidence of 0.96. 

Jackknife or bootstrap methods are often used with correlated data. These methods give esti- 
mates of the errors in fit parameters, but they do not provide information about goodness of fit. 

Once the hadron propagators are fit, we still need to perform chiral or continuum extrapolations. 
In these cases, it is also imperative to deal with the correlations among the fitted quantities that 
come from the same ensemble. With partial quenching these covariance matrices can become 
quite large, so it is essential to have enough configurations in each ensemble to be able to get a 
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FIG. 18 The p mass in units of r\ , plotted versus the squared pion mass. Since m\ oc m q , this is effectively 
a plot versus light quark mass. The octagons are from ensembles with a 0.12 fm, the squares from 
ensembles with a as 0.09 fm, and the bursts from ensembles with a m 0.06 fm. The decorated plus at the left 
is the physical p mass, with the error on this point coming from the error in r\. For reference, the upward 
arrow indicates approximately where the quark mass equals the strange quark mass. 



good estimate of the covariance matrix. 



C. Results for some light hadrons 

The pseudoscalar mesons are special for several reasons. First, very accurate mass computa- 
tions are possible. This is because the statistical error in the correlator (square root of the variance) 
decreases with the same exponential as the correlator itself - the fractional error is nearly indepen- 
dent of t, and accurate correlators can be computed out to the full extent of the lattice. Second, for 
equal mass quarks the pseudoscalar correlator does not have oscillating contributions from oppo- 
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site parity particles, and the oscillating contributions are negligible for the kaon. Third, because of 
the pions' role as the approximate Goldstone bosons for broken chiral symmetry, the breaking of 
taste symmetry leads to large mass splittings among the different taste combinations. Finally, be- 
cause it is related to the decay constant of the meson, the amplitude of the pseudoscalar correlator 
is as interesting as the mass. Because of the exact U(l) chiral symmetry of the staggered quark 
action, the axial-vector current corresponding to the Goldstone (taste pseudoscalar) pion needs no 
renormalization, so the decay constants can also be calculated to high precision. For these reasons, 
discussion of the light pseudoscalar mesons is deferred to Sec. |VH 

For the vector mesons, the fractional statistical error in the correlator increases as e^ Mv ~ Mps ^. 
Also, the vector mesons decay strongly. On the lattice, conservation of momentum and angu- 
lar momentum forbids the mixing of a zero-momentum vector meson with two zero momentum 
pseudoscalars, so the vector meson is "stable on the lattice" for pion masses large enough that 



2^JMp S + [2%/L) 2 > My. (Taste breaking adds some additional complications to this.) For all 
of the asqtad ensembles except those with the smallest quark masses, this condition is satisfied, 
and the vector meson masses can be easily, if not accurately, found. However, the problem of 
extrapolation through the decay threshold to the physical quark mass has not been fully addressed. 
Figure [18] shows the p meson mass as a function of light quark mass for three different lattice 
spacings. Results for the K* and (]) are similar, except that there is an added complication in that 
the mass needs to be adjusted to compensate for the fact that the strange q uark mass used in the 



cor relator computations differs from the physical m s . While the values in 



Bernard et al. 



and 



Aubin et al. 



(2001) 



(|2004al) use the same valence and sea strange quark masses, the masses in Fig 
have been interpolated to the correct valence strange quark mass. 

The nucleon is stable and chiral perturbation theory is available to guide the extrapolation in 
quark mass. However, computation of reliable masses is difficult because the fractional error in the 
nucleon propagator increases as e^ MN ~^ Mps ^ . Also, there are excited states with masses not too far 
above the nucleon mass that contribute to the correlator. In fact, with staggered quarks the simplest 
baryon source operators couple to the A as well as the nucle on, so the lowest positive-parity excited 
state in the correlator is the A jGolterman and Smil 1 1 9 8 sf) . Figure [191 shows nucleon masses for 
three lattice spacings versus quark mass, together wit h a continuum and chi r al extr apolation. 

Another hadron of particular interest is the Q. Ilbussaint and Davies , 2005 ). This particle 
is stable against strong decays. Also, in one-loop chiral perturbation theory there are no pion- 
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FIG. 19 The nucleon and a chiral fit. Nucleon masses are shown for different light quark masses at 
three lattice spacings. The cross at the left is the experimental value. The slightly curved line and the 
diamond at the physical quark mass are a continuum and chiral extrapolation. Lattice spacing errors are 
assumed to be linear in a 2 a y . The particular chiral form used h ere is a one-loop calculation with n — N and 



71 — A intermediate st ates (Bernard et al. 



Bernard et al. 



(I2007ch . 



1993; 



Jenkins! . 



1992). This plot is an updated version of one in 



baryon loops, so at this order there are no logarithms of m % in the chiral extrapolation of the mass. 
Therefore, we expect that a simple polynomial extrapolation in light quark mass should be good. 
Unfortunately, the Q. is a difficult mass computation with staggered quarks, first because it is a 
heavy particle and second because a bar yon operator that has t he Q~ as its lowest energ y state has 
its three quarks at different lattice sites (IGolterman and SmitL Il985l : iGupta et all Il99ll) . The Q. 
mass is strongly dependent on the strange quark mass, and in principle provides an independent 
way to determine the correct lattice strange quark mass. Figure l20l contains Q. mass estimates, 
using strange valence quark masses at each lattice spacing that were independently determined 
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FIG. 20 The Q. mass. Results are shown for three different lattice spacings. The points with a « 0.09 
fm and a w 0.06 fm were fit to the form M^r\ = A + Ba 2 a s + C(m n r\) 2 . The sloping lines show this fit 
form evaluated at the values of a 2 a s for these lattice spacings, and at a = 0. Finally, the fancy cross with 
eiTor bars is the fit form evaluated at the physical pion mass, and the small diamond is the experimental 
value. Note that in this case the ver tical axis does not beg in at zero. Earlier versions of the plot appeared in 



Toussaint and Daviesl ( 2005 ) and in 



Bernard et al. 



(I2007ch . 



from the pseudoscalar meson analysis in Sec.|VTl To do this, correlators were generated using 
two different strange quark masses near the desired one, and the f2~mass was obtained by linearly 
interpolating to the strange quark mass determined separately. This plot also shows a continuum 
and chiral extrapolation using the simple form M^r\ = A + Ba 2 a s + C{m n r\) 2 . 

Masses of othe r particles, su c h as t he a\ and b\ and particles including strange quarks 



were calculated in 
fied in 



ied in 



Bernard et al. 



Bernard et al 



Bernard et al. 



(2003b 



(2001, 



2007 d) . and the excited state of the pion was identi- 



(2007qV Light hybrid mesons with exotic quantum numbers were stud- 



cj), and exotic hybrid mesons with nonrelativistic heavy quarks in 
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Burch and Toussainti (120031) . and lBurch et g/.l (2001. 2002) 



D. Flavor singlet spectroscopy 



Determining the masses of flavor-singlet mesons is, perhaps, the most challenging endeavor in 
lattice QCD light hadron spectroscopy. The difficulty in doing so has three main sources: 

(i) Flavor-singlet correlators have two different contributions: quark-line connected and quark- 
line disconnected. The quark-line disconnected piece requires so-called "all-to-all" correlators. To 
avoid the O(V) i nversi ons to compute these all-to-all propagators, stochastic methods are used. 



Kuramashi et al 



(| 19941 ) used a unit source at eac h site and let gauge invariance do the averaging . 



I inn I ^\ 

More common now is the use of random sources (|Dong and LiuLll994l : IVenkataraman and Kilcun . 



19971) similar to Eqs. (l62l). (1631). with various noise reduction techniques (IFo 



Mathur and Dong , 



2003 



McNeile and Michaell. 



2001 



Struckmann et al. 



includ ing low-eigenmode preconditioning (IDeGrand and Hellei 
19981) . 



2002 



2001 



ev et al 



Wilcox, 



2005 



19991) 



Venkataraman and Kilcup 



(ii) While the stochastic noise of the quark-line connected correlators falls off exponentially 
(albeit with a smaller exponent than the signal), the noise in the quark-line disconnected part is 
constant. So the signal to noise ratio falls off much faster for the disconnected part. 

(iii) The quark-line connected correlator is the same as for a flavor-nonsinglet meson - in partic- 
ular the pion for the pseudoscalar channel. Therefore, the very noisy disconnected correlator first 
has to cancel the connected correlator before giving the desired singlet correlator whose falloff 
gives the flavor-singlet mass. 

Since much larger statistics are needed for the computation of the flavor- singlet correlators, 



the UKQC D collaboration 



trajectories ([Gregory et al. 



las extended a co uple of the MILC lattice ensembles to around 30000 



2001 



2008, 



20091) . Their simulations are still on-going. So far, the only 



result given is for the ++ glueball, whose correlator can be constructed from gauge field operators 
and requires no noisy estimators and Dirac operator inversions. For two different lattice spacings, 
a ^0.12 and . 09 fm , the UKQCD collaboration finds m 0++ = 1629(32) MeV and 1600(71) MeV 



(IGregory et al 



20091) . respectively. 



It is important to continue this investigation. In particular, obtaining the correct x\ r mass would 
further support the correctness of the rooting procedure to eliminate the unwanted tastes for stag- 
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gered fermions. 



E. Scalar mesons fo and a 



In this subsection, we describe briefly the analysis of correlators for two light, unstable scalar 
mesons, namely, the isosinglet fo and the isovector ao- 

With the first good measurements of the ao channel in the staggered fermion formulation a 
peculiarity was encountered: it was found that on coarse lattices the ao correlator appeared to have 
a spectral contribution with an anomalou sly low mass, lighter than any physical decay channel 



(Aubin et al. 



2004a 



Gregory et al. . 



2006) 



For sufficiently light u and d quark masses, the fo decays to two pions. Likewise, the isovector 
scalar meson ao decays to a pion and an T|. On the lattice, the open decay channels complicate 
the analysis of the scalar meson correlators. They are dominated by the spectral contributions of 
the significantly lighter decay channels. As a flavor singlet, the fo also suffers from the quark-line 
disconnected contributions described in the previous subsection. Finally, with staggered fermions 
at nonzero lattice spacing, the splitting of the pseudoscalar meson taste multiplets in the decay 
channel deals a seeming coup de grace. 



Fortunately, on e can make progress using rSXPT described in Sec. lIII.AkBernard et all 



Prelovsekl 



2006a; 



2006allb ). The essential idea is to match definitions of the desired correlator of local 



interpolating operators in the lattice QCD formulation and in rSXPT. The lattice definition is the 
basis for the numerical simulation of the correlator, and the rSXPT definition provides a model for 
fitting the result of the simulation, including all taste-breaking effects in the decay channels. If we 
take the taste-multiplet masses from separate calculations, then, despite the rather complicated set 
of two-meson channels, that portion of the fit model depends on only three low energy constants. 
In principle, even these constants can be determined from independent measurements, leaving no 
free parameters. So this fit provides a further test of the viability of rSZPT as a low energy effective 
theory for the staggered action. 

The hadron propagator from lattice site to y is defined in the same way from the generating 
functionals for both QCD and the chiral theory: 

a 2 logZ 



dm f j>(y)dm e/ , e (0) 



(128) 
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FIG. 21 Best fit to the ciq correlator (left panel) for five momenta and the fo correlator (right panel) for four 
momenta. The fitting range is indicated by points and fitted lines in red and blue (darker points and lines). 
Occasional points with negative central values are not plotted. Da ta are determined from the a»0.12 fm 



(coarse) ensemble with ami = 0.005 and am s = 0.05. Figures from 



Bernard et al. 



(12007 ah . 



In QCD, the source mt f(y) generalizes the usual quark mass term and includes off-diagonal flavor 
mixing /,/'. The same correlator is defined in rSZPT, where the local source rnfj'iy) appears 
in the generalized meson mass matrix. This establishes a correspondence between the correlator 
defined in terms of the quark fields q(y)q(y) in QCD and in terms of the local meson fields B<t> 2 (y). 

To lowest order in rSXPT, the meson correlator is described by a bubble diagram, which gives 
the contributions of the two-pseudoscalar-meson intermediate states, including all taste multiplets 
and hairpins. These contributions are determined from the multiplet masses and the rSZPT low 
energy constants B, d' A , and d' v described in Sec. IIII.Al In addition to the bubble diagram, one adds 
an explicit quark-antiquark ciq or /o state to complete the fit model. Results are shown in Fig.l2Tl 
and results for the low energy constants are listed in Table ITVl 

It is particularly instructive to examine the variety of two-pseudoscalar-meson taste channels 
contributing to the scalar meson correlators. To be physical states, the external scalar mesons ao 
and fo must be taste singlets. Taste selection rules then require that they couple only to pairs of 
pseudoscalar mesons of the same taste. Thus, for example, for the ao, each flavor channel, such 
as 71 — T|, comes with a multiplicity of sixteen taste pairs, although lattice symmetries reduce the 
number of distinct thresholds to six. There is also a set of n — r\' channels. To get the energies 
of the thresholds, we look at the taste splitting of the component hadrons. We have already seen 
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fo and ao correlators Meson masses and decays 


riml/(2m Utd ) 
5a 


7.3(1.6) 6.7 
(prior) -0.016(23) 
-0.056(10) -0.040(6) 



TABLE IV Comp arison of our fit parameters for the rSXPT low energy constants with results from 



Aubin et al. 



(2004b) 



how the pion taste multiplet splits into the Goldstone state and a variety of higher- lying states, 
all of which become degenerate in the continuum limit. The T| and r\', on the other hand, have 
unusual splitting because they mix with the chiral anomaly. Since the anomaly is a taste singlet, 
only the taste-singlet T| and r\ f mix with it in the usual way. Thus, in the continuum limit only 
the taste singlet states are expected to have the correct masses. They are the only physical states. 
The fifteen taste nonsinglet T|'s and r^'s remain light. The pseudoscalar-taste eta pairs with the 
pseudoscalar-taste pion. The unphysical pse udoscalar-taste 7t — r\ channel gives an anomalously 
light spectral contribution to the ao correlator (IPrelovsekll2006alJbh . A similar complication occurs 
in the fo correlator, but it is masked by the expected physical two-pion intermediate state. 

The unphysical taste contributions provide a concrete illustration of the breakdown of unitarity 
at nonzero lattice spacing as a result of the fourth-root. The theory heals the scalar meson cor- 
relators in the continuum limit by a mechanism that parallels exactly the one described for the 
one-flavor model in Sec. IIII.Cl The pseudoscalar meson bubble diagram contains a negative-norm 
channel. This unphysical ghost channel has the weight needed to cancel the contributions of all the 
unphysical taste components in the continuum limit. Thus in the continuum limit only the physical 
intermediate two-meson states survive. 

The behavior of the isovector scalar correlator has al so been analyzed fo r the case of domain- 



2008). In the mixed-action 



wall valence quarks on the MILC staggered ensembles (lAubin et al. , 
case, the ao correlator receives contributions from two-particle intermediate states with mesons 
composed of two domain-wall quarks, mixed mesons composed of one domain- wall and one stag- 
gered quark, and mesons composed of two staggered quarks. Because the symmetry of the external 
valence quarks restricts the sea-sea mesons to be taste singlets, the correlator does not receive con- 
tributions from all of the taste channels. As in the purely staggered case, the one-loop bubble 
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FIG. 22 The isovector scalar (ao) correlator on the MILC coarse ami/am s = 0.007/0.05 ensemble with 
three different domain-wall valence masses. Overlaid on the data are the predicted bubble contributions, 
whic h should dominate over the exponentially-decaying contributions at sufficiently large times. Figure 



from 



Aubin et al. 



(2008). 



2006bl) . which are known 



contribution is determined by three low-energy constants (IPrelovsekl 
from tree-level XPT fits to meson masses. For domain- wall quarks on the coarse and fine MILC 
lattices, the contribution from the bubble term is predicted to be large and negative for several time 
slices. Thus a comparison of the mixed-action XPT prediction for the behavior of the ao correlator 
wi th lattice data provi des a strong consistency check. 



Aubin et al. 



(|2008l) compare the mixed-action XPT prediction for the bubble contribution with 
the lattice ao correlator for several domain- wall valence masses on the coarse and fine MILC 
lattices. They find that, in all cases the size of the bubble contribution is quantitatively consistent 
with the data, and that the behavior of the data cannot be explained if mixed-action lattice artifacts 
are neglected. For fixed light sea quark mass, the size of the bubble term decreases as the valence 
quark ma ss increases (see F ig. 1221) . The bubble contribution also decreases as a — ► 0. These 



results of 



Aubin et al. 



(|2008r) support the claim that mixed-action XPT is indeed the low-energy 



effective theory of the domain-wall valence, staggered sea lattice theory. Furthermore, mixed- 
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action %PT describes the dominant unitarity-violating effects in the mixed-action theory even when 
such effects are larger than the continuum full QCD contributions that one wishes to extract. Thus 
mixed-action %PT fits can be used to remove taste-breaking and unitarity-violating artifacts and 
recover physical quantities. 



F. Summary 



In general these and other lattice spectrum calculations confirm that QCD does predict the 
hadron spectrum. However, although we can see the effects of decay thresholds as the quark mass 
is varied (e.g., Sec. IV.EI) . a nd though some s cattering lengths can be indirectly determined through 



chiral perturbation theory (|LeutwylerL 
to be calculated in the future. 



2006f) . most hadronic decay rates and cross sections remain 



VI. RESULTS FOR THE LIGHT PSEUDOSCALAR MESONS 



A. Motivation 



Precise computations are possible for light pseudoscalar mesons (see Sec. IV.CI) . and they lead 
to interesting physics. If lattice calculations of light pseudoscalar mesons and decay constants 
can approach the chiral and continuum limits, we can determine the up, down and strange quark 
masses and many of the low energy constants (LECs) of the chiral Lagrangian, includ ing several 



19841) . From the 



combinations of the NLO Gasser-Leutwyler constants L, (IGasser and Leutwylerl . 
ratio /k/I-k, we can extract \V US \ from the kaon leptonic branching fraction, providing a test of 
CKM matrix unitarity for the first row of the matrix. 



B. From correlators to lattice masses and decay constants 



Study of the light pseudoscalar mesons on MILC lattices beg an in 2004 (lAubin et al 



and has included several updates at th e annual Latt i ce conf erences ([Bernard et al. 



We first review the methodology of 



Aubin et al. 



2006d. 



2004b) 



2007eh 



(|2004br) . In the Goldstone (taste pseudoscalar) 



case, we can use the PCAC relation to relate the decay constant fp$ to matrix elements of the spin- 
and taste-pseudoscalar operator Op(t) = \j/(y5 0^5 )\|/ between the vacuum and the meson. In terms 
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of the one-component staggered quark formalism, 

P (t) = % a (x,t)(-lf + Y(x,t), (129) 
where a is the (summed) color index. As in Eqs. (112111 22b . we define a correlator by 

Cpp(t) = ^Y,(°P(y,t)ol(x,0))=c P pe- m ^ + ... , (130) 



y 



where mps is the mass of the (lightest) pseudoscalar and V s is the spatial volume. After fitting the 
correlator to this form, we can find the decay constant from 



^ = K + m y )^/^g, (131) 

where m x and m y are the two valence quark masses. 

Although the decay constant is found from the overlap of the point-source operator with the 
meson state, most directly obtained from the point-point correlator Eq. (11301) . it is useful to use the 
Coulomb wall source Eq. (ll24l) and point sink to calculate the correlator 

C W p = (Op(x,t)ol(0))=c WP e- mpst + ... . (132) 

The advantage of this correlator is that it has less contamination from excited states than does Cpp, 
and helps in fixing the pseudoscalar mass. 

A random-wall source can also be used instead of a point source to calculate Cpp, giving smaller 
statistical errors. The source for the quark on each site of a time slice is a three component complex 
unit vector with a random direction in color space. Thus, contributions where the quark and anti- 
quark in a meson originate on different spatial sites average to zero. After dividing by the spatial 
lattice volume, this source is used instead of Op in Cpp. The preferred method is then to fit Cwp 
and the random-wall point-sink Cpp with three free parameters App, A\yp and mps'. 

Cpp = m\ s A PP e- m ™\ 

C WP = ml s A WP e- mps \ (133) 

so that App is the desired combination cpp/mp s that appears in Eq. (11311) . An appropriate range of 
Euclidean time must be selected to get a good confidence level of the fit. If the minimum distance 
from the source point is too small, there will be excited state contamination. It is essential to use 
the full correlation matrix of the data to get a meaningful confidence level and avoid contamination. 
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For chiral fits used to extract LECs that govern the mass-dependence of physical quantities, it 
is important to fix the scale in a mass -independent manner. This is because all mass dependence 
should be explicit in %PT, and none should be hidden in the scale-fixing scheme. As described in 
Sec. IIV.BI a mass independent method is used to determine a in which r\/a is extrapolated to the 
physical, rather than simulated, quark masses on the given ensemble. 

Partial quenching is very useful in order to obtain enough data to perform the required chiral 
fits. For the valence masses on a typical ensemble, nine different masses from 0.lm' s to m' s (m' s is 
the simulated strange sea mass) may be used. This yields 45 distinct pairs of valence masses, and 
hence 90 values (meson masses and decay constants) for the chiral fit. Without partial quenching, 
we would have only four values. Of course, the correlations among the 90 values must be taken 
into account. 

Finite volume corrections are included in the one-loop rSZPT forms used to fit the lattice data. 
Since the spatial box sizes are at least 2.4 fm, and for the smallest light sea-quark masses they are 
increased to about 2.9 fm or larger, these corrections are always less than 1 .5%. Smaller, additional 
corrections representing "residual" effects from higher-loop contributions are applied at the end of 
the calculation, as described below. The results cannot be fit without the on e-loop finite volume 



Aubin et al, 



(|2004bh . five coarse and two 



corrections, nor can they be fit with continuum %PT. In 
fine ensembles were fit with continuum XPT; however, the confidence le vel of the fit was 10~ 2 50 ! 
In the remainder of this section, we present methods and results from 



Bernard et al. 



(I2007eh. A 



final version of the analysis, using add ed ensembles and two-loop chiral logarithms (|Bijnens et all 



2004, 



2006; 



Bijnens and Lahdel . 120051) . is in progress. 



The fitting is done in two stages. In the first stage, the leading order (LO) and next-to-leading 
order (NLO) low energy constants (LECs) are determined by fitting a restricted set of data that 
is closer to the chiral and continuum limits than the additional points included later. Specifically, 
the largest lattice spacing (a ~ 0.15 fm) is omitted and the valence quark masses are required 
to obey am x + am y < 0.39 am s (for a ~ 0.12 fm), am x + am y < 0.51 am s (for a ~ 0.09 fm), and 
am x + am y < 0.56 am s (for a ~ 0.06 fm). Further, for a ~ 0. 12 fm three higher- mass combinations 
of sea-quark masses are omitted. Despite the restrictions, it is found that due to the high precision 
of the data it is necessary to add NNLO analytic terms in order to get good fits. In the second 
stage, the range of valence and sea-quark masses is extended to include the region around the 
strange quark mass. The LO and NLO low energy constants are constrained to be within the range 
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FIG. 23 NNNLO fit to partially-quenched squared meson masses. Only the lightest sea-quark ensemble for 
each lattice spacing is shown. The data fit includes the results fo r decay constants and is reflected in the 



number of degrees of freedom. Figure from 



Bernard et al. 



(I2007eh . 



determined by the first stage of fitting. In this stage, NNNLO analytic terms are needed to get good 
fits. 

In Fig.[23l we show the squared meson masses in units of (GeV) 2 . For the "pions" m x = m y . For 
the "kaons" a few fixed values of m y are picked for illustration, and m x is varied. The horizontal 
axis is m x /m' s . Only a small fraction of the points used in the fit are shown. For each lattice spacing, 
the plot contains only the lightest sea-quark mass ensemble, and no decay constant data is plotted. 
For this fit, % 2 = 436 with 449 degrees of freedom, corresponding to a confidence level of 0.66. The 
dashed red line shows the continuum prediction after all lattice spacing dependence in the fit pa- 
rameters is extrapolated away, the strange sea-quark mass is fixed to its physical value and the light 
valence and sea masses are set equal. The physical values of m s and m = (m u + ma)/2 are required 
to simultaneously yield the kaon and pion masses denoted K and ft in the figure. These masses 
correspond to what the kaon and pion masses would be with isospin and electromagnetic effects 
removed. Some pheno menological input is n eeded to account for the electromagnetic effects. This 



is explained in detail in 



Aubin et al. 



(|2004bl) . The vertical do tted line is drawn at m/m 



The "residual" finite volume corrections are then applied. IColangelo. Purr, and Haefelil (120051) 
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FIG. 24 The meson decay constants are plotted along with the NNNLO fit that was shown for the masses 
in Fig. [23] The left plot shows partially quenched data from more ensembles than in Fig. |23l but still only 
a fraction of the data fit. On t he right, still more ens embles are included, but only full QCD data points are 



plotted. Both figures are from 



Bernard et al. 



(12007eh . 



have shown that higher than one-loop %PT corrections can be significant in the current range of 
quark masses and volumes. For a w 0.12 fm with sea masses ami/am' s = 0.01/0.05, there is a 



direct test of finite v o 
sides. Bernard et al. 



ume ef fects on 20 3 and 28 3 volumes that correspond to 2.4 and 3.4 fm box 



(|2007e[) detail the comparison between these calculations and the one-loop 



result. On this basis, a small correction is applied to the continuum prediction. This amounts to 
0.25% for /„, 0.05% for f K , -0.15% for n%, and -0.10% for m\. These values are also added to 
the systematic error. 

By extending the kaon extrapolation line in Fig.|23j one finds the value of m u that corresponds 
to the K + mass (see 



Aubin et al. 



(|2004b|) ). Two important mass ratios are determined: 



m s /m = 27.2(1)(3)(0) , m u /m d = 0.42(0)(1)(4) . 



(134) 



The errors are statistical, lattice-systematic, and electromagnetic (from continuum estimates). Note 
that the m u = solution to the strong CP problem is ruled out at the 10 a level. 

Having found the continuum fit parameters and the physical quark masses, the decay constants 
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are predicted. Figure [24] (left) shows (some of) the decay constant data and the fit through the 
displayed data. For the continuum prediction (dashed red line), the strange sea-quark mass is set 
to its physical value and the light valence and sea masses are set equal. The left end of the curve 
corresponds to m x = m y — m. The vertical error bar to the left of the + shows the systematic error. 



The experimental result is shown as an octagon. 
assumption that |V^| = 0.97377(27) (|Amsler et al. 



t come s from the decay n + — > /j + V jU with the 



20081) . Figure [21 (right) shows the full QCD 



points from a slightly different fit with data from additional ensembles. Note that the data points 
at a ~ 0.06 fm are quite close to the full QCD continuum extrapolated curve. 

Up to this point, the lattice spacing is set by calculation of the heavy-quark potential parameter 
r\, which yields relative lattice spacings between ens embles, and the co ntinuum extrapolation of 



T splittings determined by the HPQCD collaboration (jGray et al. 
scale. These results yield a value r\ = 0.318(7) fm. On this basis, 



2005|), which gives an absolute 



fn 
h 
In/In 



128.3 ±0.5 1^ MeV : 
154.3 ± 0.4 tH MeV, 
1.202(3)(+ J) , 



(135) 



where the errors are from statistics and lattice system atics. This value for f n is consistent with the 



2008) 



experimental result, f^ xpt = 130.7 ± 0. 1 ± 0.36 MeV dAmsler et all 

An alternative approach is to set the scale from f % itself. In this case, there are small changes in 
the quark masses and 

n =0.3 108 (15) (^9) fm, (136) 

which is 1-G lower (and with somewhat smaller errors) than the value from the T system. For the 
decay constants, 



f K = 156.5 ± 0.4 t\j MeV : 
/*//* = U97(3)CiS 



(137) 



w here the errors are statistical and systematic. 



Marciand (|2004l) noted t hat the lattice value of fx] fn can be combined with measurements of 



the kaon branching fraction (|Ambrosino et all l2006aUbl) to obtain |V ; „|. From Eq. (11371) 



|V MS |=0.2246(+?f) 



(138) 
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which is consistent with (and competitive with) the world- average value \V US \ = 0.2255(19) 



(Amsler et al. 



2008) coming from semileptonic A"-decay coupled with non-lattice theory. 
Using the two- loop perturbative calculation of the mass renormalization constant Z„ 



(|Mason et al. 



2006I) 11 . absolute quark masses can be found: 



m s = 88(0)(3)(4)(0) MeV 
m M = 1.9(0)(l)(l)(l)MeV 



m = 3.2(0)(l)(2)(0) MeV, 
m rf = 4.6(0)(2)(2)(l)MeV. 



(139) 



The errors are statistical, lattice-systematic, perturbative, and electromagnetic (from continuum 
estimates). Nonperturbative computations of Z m are in progress. 

The chiral fits also determine various Gasser-Leutwyler low energy constants and chiral con- 
densates: 



2L 6 — L\ 


= o.4(i)(; 


1), 


2L 8 -L 5 = -0.1(1)(1), 


U 


= 0.4(3)(; 


1), 


Z* = 2.2(2X1?) , 


U 


= 0.4(2)(; 


1), 


L 8 = 1.0(l)(l), 


hi Si = 


1.052(2)(; 




( a - M ) 2 = -(278(l)(^)( 5 )MeV) 3 . 


fn/fs = 


1.21(5X1 




( M - M ) 3 = -(242(9)(l 1 5 7 )(4) MeV) 3 


h/h = 


1.15(5X1 




(uu) 2 / (m) 3 = 1.52(17)(±f 5 ) . 



(140) 

The errors are statistical, lattice-systematic and perturbative for the condensates. f 2 (^3) denotes 
the three-flavor decay constant in the two (three) flavor chiral limit, and (uu)2 ((uu) 3 ) is the corre- 
sponding condensate. The low energy constants L,- are in units of 10~ 3 and are evaluated at chiral 
scale win; the condensates and masses are in the MS scheme at scale 2GeV. The indications are 
that the L,- will change signifi cantly when th e two-loop logarithms are included, just as they do in 



phenomenological estimates (|Bijnens 



20071) . Other results are very stable, however. 
The rSXPT formalism relies on the replica proceedure, and taking the fourth root corresponds 
to n r = 1/4 where n r is the number of replicas. The fact that there are good fits with the rSXPT 
formulae, but not with continuum XPT, is a test of staggered chiral perturbation theory. A further 



11 With this two-loop Z,„-factor a tadpole improved definition of the bare quark mass should be used, in which what 
we have denoted by am q throughout this review should be replaced by UQ<mi q . The tadpole factors for the various 
MILC ensembles are listed in Table U 
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test of rSXPT is to allow n r to be a free parameter in the fits. For the low mass data, n r = 0.28(2) (3) 
where the first error is statistical and the second systematic coming from varying the details of 
the chiral fits. We are encouraged by this strong constraint on n r , and the success of rSXPT in 
describing the MILC data. 



C. Other computations of f n and fx 



Since the MILC collaboration's initial calculation of the light pseudoscal ar meson masses, de 



2004b). 



cay constants, and quark masses using the a ~ 0.12 and 0.09 fm lattices (|Aubin et al. . 
several other groups have also computed f n and /k on the MILC ensembles using different va- 
lence quark formulations. All of the results are consistent with those of the MILC collaboration, 
Eq. (11351) . and with each other. 

The HPQCD collaboration uses HISQ staggered valence quarks and the MIL C asqtad stag- 
gered sea quark ensembles with lattice spacings a « 0.15, 0.12 and 0.09 fm fm (IFollana et all 
20081) . They generate one "pion" point and one "kaon" point per ensemble, matching the masses 
of the Goldstone HISQ pi on to the asqtad one , and the mass of the HISQ ss meson to 696 MeV, 



the XPT value. Although 



Follana et al. 



(|2008l) are performing a mixed action lattice simulation, 
they extrapolate to the physical light quark masses and the continuum using continuum NLO XPT 
augmented by analytic terms constrained with Bayesian priors. Terms proportional to a s a 2 and 
a 4 are included to test for conventional discretization errors, while those proportional to 0? s a 2 , 
<y? s a 2 \o%(m q ) , and a? s a 2 m q are intended to test for residual taste-changing interactions with the 
HISQ valence quarks. HPQCD obtains the following results for f % , fx, and the ratio: 



fn = 132(2)MeV, f K = 157(2)MeV, f K /f n = 1.189(7), 



(141) 



where the largest source of error is the uncertainty in the scale r\ (1.4% for f % and 1.1% for fx). 

The NPLQCD collabor ation uses domain-w all valence quarks and four a « 0.12 fm ensembles 
with mi/m' s = 0.14 - 0.6 (IBeane et ail 12007 ah . They tune to match the valence pion and kaon 
to the corresponding asqtad particles. Due to the mixed action, there are still unitarity-violating 
artifacts that vanish only in the limit a — > 0. They compute only the ratio fx/fn, which has a 
milder dependence upon the quark mass than the individual decay constants, and extrapolate to 
the physical light quark masses using the NLO continuum XPT expression, which depends only on 
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one free parameter, L5. The result is 



fit/ f% 



1.218±0.002±g;gJi , 



L 5 {m^) = 2.22 ±0.02+g;^ x 1(T 3 



(142) 



where the first error is statistical and the second error is the sum of systematic errors added in 
quadrature. The dominant source of uncertainty is from the truncation of the %PT expression 
(-0 022 f° r me rat i°)> which they estimate by varying the fit function through the addition of NNLO 
analytic terms and double logarithms. Although they do not include an error due to their use of 



on 



V a single la t tice spa cing, this is likely a small effect in the ratio fn/fn- 



Aubin et al. 



(I2009al) also use domain-wall valence quarks. In contrast with NPLQCD, however, 
they compute many partially quenched points on the a ~ 0. 12 and 0.09 fm ensembles, and use NLO 
mixed actio n XPT with highe r-order analytic terms to extrapolate to physical quark masses and the 



continuum (|Bar et all 
constants are 



2005). Their preliminary results for the light pseudoscalar meson decay 



A = 129.1(1.9)(4.0)MeV, f K = 153.9(1.7)(4.4)MeV, f K /h = 1 - 191 (16) ( 17) , (143) 

where the first error is statistical and the second is the sum of systematic errors added in quadrature. 
The dominant source of error is from the chiral extrapolation procedure (2.2% for f n and 2.3% for 
fx), and is estimated by varying the analytic terms included in the fit function. 

In Fig.[25l we compare results for fx;/ fn from a variety of 2+1 flavor calculations. The top four 
results all use MILC aqstad configurations and were discussed above. The two lower results from 



the PACS-CS Collaboration (Aoki et al. 



200% and the RBC/UKQCD Collaboration (lAllton et al. . 



2008|) use clover quarks with the Iwasaki gauge action and domain-wall quarks, respectively. 



VII. HEAVY-LIGHT MESONS: MASSES AND DECAY CONSTANTS 

Calculations of B- and D-meson masses and decay constants using the 2+1 flavor MILC con- 
figurations have been performed in joint work by the Fermilab Lattice and MILC collaborations, 
and by the HPQCD collaboration. Of numerical quantities involving heavy b and c quarks, me- 
son masses and decay constants are among the simplest quantities to compute numerically and 
are often well measured experimentally. Thus they provide valuable cross-checks of lattice QCD 
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FIG. 25 The ratio of light decay constants fy/fn from six calculations. The top four use MILC asqtad 
configurations and the lower two use other types of quarks. Details and references can be found in the text. 

methods. In particular, once the treatment of the light sea and valence quarks has been validated 
within the light pseudoscalar sector, calculations of heavy-light meson masses and decay constants 
allow tests of the various lattice QCD formalisms used for heavy quarks. In this section, we de- 
scribe the 2+1 flavor calculations by Fermilab/MILC and HPQCD of heavy-light meson masses 
and decay constants, and show that, with one exception, they are consistent with experiment. These 
results give confidence in other lattice QCD calculations involving b and c quarks, such as those 
of semileptonic form factors described in Sec. I Villi 

A. Heavy quarks on the lattice 

Heavy quarks, i.e., those for which the quark mass in lattice units am is large, present spe- 
cial challenges. As long as am <C 1, heavy quarks on the lattice can be treated with light-quark 
formalisms such as staggered fermions. At the lattice spacings in common use, we have am c ~ 
0.5-1.0 and ami, ~ 2-3. For charm quarks, light-quark methods can only be used if they are 
highly improved to remove discretization errors. Bottom quarks still require special heavy-quark 
methods. 



MILC 
NPLQCD 
HPQCD 
ALV 

(preliminary) 

PACS-CS 
RBC/UKQCD 
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1. Nonrelativistic QCD 



A straightforward way of formulating heavy quarks on the lattice is to rewrite the 



Dirac-like 



in HQET (llsgur anc 



ight-quark acti o n as a sum in a no nrelativistic operator expansion, as is done 



Wise, 



dCaswell and LepageLll986l : 
Snrqcd = Y,^( x . 

x 

where 



1992; 



Neubert, 



Lepage et al. 



1994b and in nonrelativistic expansions in QED 



1992) 



vJ, + VW = - {U„(x)y\f(x + afi)-y(x)) 



(144) 



(145) 



and where the (|) are two-component fermions representing the quarks. An analogous term in the 
action governs the antiquarks. The leading heavy-quark mass dependence is absorbed into the 
fermion field and vanishes from explicit calculations. For b quarks in particles with a s i ngle h eavy 



1990). In 



quark, the first term in this action yields the static approximation (|Eichten and HillL 
heavy-light systems, the importance of operators in this expansion is ordered according to HQET 
power counting (k ~ A/mg). In quarkonium systems, operators are ordered by heavy-quark ve- 
locity. 



2. Wilson fermions with the Fermilab interpretation 

In NRQCD, the kinetic energy operator of the Dirac action, \j7(x) L/Y; V ; \)/(jc) is replaced by the 
leading kinetic energy operator ^ (x) J- J^,-A, (|)(x) plus a series of higher dimension operators. The 
action for Wilson fermions contains the leading kinetic energy operators of both the Dirac and the 
nonrelativistic actions, as in Eq. (031) : 

x \ M V ) 

The effects of the Laplacian term, which eliminates the doubler states, vanish in the limit am — > 0. 
As am becomes larger, the importance of the Laplacian term grows. When am ^> 1, the Lapla- 
cian term dominates the Dirac-like kinetic energy term, and the theory behaves like a type of 
nonrelativistic theory in which the rest mass mi = E(p 2 = 0) does not equal the kinetic mass 
m2 = \/{2dE/dp 2 ). (Note that we use lower-case m to refer to quarks and capital M to refer to 
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mesons in this section.) As am — > 0, the two masses converge to the bare quark mass m. For 
heavy quarks the kinetic mass controls the physics, and the rest mass may be absorbed into a field 
redefinition. This means that the Wilson action and related actions can be used as actions for 
heavy quarks as long as m2, with contributions from bot h terms in the kinetic energy, is adjusted to 
equal the desired physical mass (|E1-Khadra et a/1 1 19971) . It is possible to set m\ = m.2 by breaking 
time-space axis-interchange symmetry in the Lagrangian. If this is not done, m\ and /«2 have the 
tree-level form 

am\ = log(l +amo) (147) 

and 



1 



1 



7o~Z S + TT • (148) 

am2 amo (2 + am,Q) 1 + artiQ 

The action of the nonrelativistic expansion can be viewed as arising from a field transformation 
of the Dirac field, the Foldy-Wouthuysen-Tani (FWT) transformation. The Wilson action, with 
both types of kinetic energy operators, can be viewed as arising from a partial FWT transforma- 
tion. Like the action of NRQCD, it produces the same physic s as the Dirac action as lon g as a 



2008). The 



series of correction operators is added to sufficient precision (|Oktay and KronfeldL 
leading dimension-five correction operator has the same form for heavy Wilson fermions as for 
light cloverAVilson fermions [Eq. (fT9ll. Ssw = ^^sw £ x \f/(x)G^ v ^v(*)y(*)- All simulations to 
date using this approach to heavy quarks have therefore used cloverAVilson fermions. A systematic 
improvement program is possible as outlined in Sec. IX.CI 



3. The HISQ action 



Because 0.5 <am c < 1 at currently accessible lattice spacings, it is possible to use ordinary 
light-quark actions to treat the charm quark. However, to obtain high precision it is necessary to 
correct t he action to a high or der in am. This approach is followed with "highly improved staggered 
quarks" (|Follana et q/ll2007|) . as explained in Sec. III.EI 
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B. Lattice calculations of masses and decay constants 



As in the light pseudoscalar meson case, the heavy-light decay constant is proportional to the 
matrix element of the axial current: 



(0\A t ,\H q (p)) = if H n fl , 



(149) 



where = qj^JsQ. Because of the heavy-quark normalization in HQET, it is often useful to 
consider the combination decay amplitude 



§H q =f Hq JM Hq 



which is computed from the correlators 



Cb(f) = (o Hq (t)oUo)), c M {t) = (A 4 (t)oU0)). 



(150) 



(151) 



For the case of Fermilab heavy quarks or NRQCD b quarks, the heavy-light meson mass is obtained 
from the kinetic mass (A/2) in the dispersion relation, whereas for HISQ charm quarks Mi = A/2, 
so both are simultaneously set to the D- or ZV meson mass. 

The Fermilab Lattice and MILC collaboration s' calculation of heavy-light meson decay con- 



stants (Aubin et al. 



2005a 



Bernard et al. 



2009bl) employs the Fermilab action for the heavy b 



and c quarks and the asqtad staggered action for the light u, d, and s quarks. They construct the 
heavy-light meson interpolating operator and axial vector current A M using the meth od for combin 



ing fou r-component Wilson quarks with 1 -component stagg ered quarks described in lWingate et al. 

2009bl) uses data on the 



(120031) . Their most recent determination from Lattice 2008 (IBernard et all 
medium-coarse, coarse, and fine lattices, with 8-12 partially quenched valence masses per ensem- 
ble. The clover coefficient csw and hopping parameter K in the Fermilab action are tuned to remove 
errors of O (1/mg) in the heavy-quark action. In particular, they set csw = Uq 3 , the value given by 
tree-level tadpole-improved perturbation theory (|Lepage and Mackenzie!. 1 19931 ). They choose the 
charm quark hopping parameter K c so that the spin-averaged (kinetic) D r meson mass is equal to 
its physical value, and choose the bottom quark hopping parameter K^, to reproduce the 5 v -nieson 
mass in an analogous manner; this implicitly fixes the b and c quark masses. They also remove 
errors of C(l/mg) from the heavy-light axial vector current A^ by rotating the heavy-quark field 



in the two-point correlation function: 



Mf b -»■ m b = 1 +ada- D M \\f b 



(152) 
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where D\ ai is the symmetric, nearest-neigh 



?or, covariant d i fferen ce operator, and the tadpole- 



improved tree-level value for d\ is given by (|E1-Khadra et al. 



1997): 



d\ = — 



1 



+ 



1 



(153) 



wo \2 + amo 2(l+amo) 
They obtain the renormalizatio n factor needed to mat ch the lattice heavy-light current onto the 



continuum using the method of 



Hashimoto et al, 



(1999) 



z9 q = p% q Jz£ Q z% q 



-A 4 " (154) 
where the flavor-conserving factors Zy® and Z y q are determined nonpe rturbatively and t he re 



maining factor is determined to 1-loop in lattice perturbation theory (lEl-Khadra et al. 



Lepage and Mackenzie . 



1993) 



2007; 



The Fermilab/MILC collaboration fits its decay constant data as a function of light-quark sea 
and valence masses to the one-loop form given by HMSXPT (see Sec. IIII.BI) . supplemented by 
analytic NNLO terms, which are quadratic in the light valence and/or sea masses. This is very 
similar to the approach taken in the light pseudoscalar sector, as described in Sec.|VlJ While pure 
NLO fits are adequate to describe the data for very light valence mass, once this mass gets to be 
roughly half the strange quark mass or higher, at least some NNLO terms are necessary to obtain 
acceptable fits. 

Figure [26] shows the preferred HMS5CPT fit to data at multiple lattice spacings for 4>o and 4>d s , 
which are functions of the light valence mass, the light sea mass and the strange sea mass. In 
addition to taste-breaking discretization effects that appear as taste-splittings, taste-hairpins, and 
taste-violating analytic terms, there are "generic" light-quark discretization effects, which can be 
thought of as changes in the physical LECs (such as 4>o, the value of 4> in the SU(3) chiral limit) 
with lattice spacing. With the asqtad action, such effects are O (asa 2 ). They can be (approximately) 
accounted for by adding additi onal paramete r s to th e HMSXPT fit function, with variations limited 
by Bayesian priors, following [Lepage et all (|2002l) . This is done in the fit shown in Fig. [26l al- 
though the effects are quite small, and fits without the additional parameters give almost the same 
results (and confidence levels), but with somewhat smaller statistical errors. Once the parameters 
of the HMSXPT fit are known, taste- violating and generic discretization effects through O (a 2 ) can 
be removed by setting a = 0. After taking the continuum limit, the valence and sea-quark masses 
are set to their physical values in order to obtain the decay constants of a D + and D s meson, up to 
tiny isospin violations in the sea sector. 
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FIG. 26 Chir al extrapolation for 



collaboration (IBernard et al. 



<3>d (octagons) and <J>o s (crosses or diamonds) by the Fermilab/MILC 
2009bl) . Solid lines are the HMSXPT fit to <& D ; dotted lines, to <Z> Ds . Although 
the full set of partially-quenched data is included in the fit, for <t>o the plot shows only those full QCD points 
for which the light valence and sea masses are equal to m x , the mass on the abscissa. For <£>d s , only points 
with the strange valence mass (m sv ) equal to the strange sea mass are shown, plotted either as a function 
of mi (crosses), or at m sv (diamonds). The (red) dashed lines show the fit after removal of light-quark 
discretization errors, with the fancy plus signs giving the chirally extrapolated results with statistical errors. 



The HPQCD collaboration's calculation of the B and ^-meson decay constants ( 



Gamiz et ail 



2009|) employs the NRQCD action for the heavy b quarks and the asqtad staggered action for the 
light u, d, and s quarks. They use six data points in their analysis — four full QCD points on the 
coarse ensembles and two full QCD points on the fine ensembles. They fix the Z?-quark mass so 



that the mass of a bb meson reproduces the physical my (IGrav et al. 



tion includes all currents of 0(1/ mi,) (|Morningstar and ShigemitsuL 



perturbation theory to match onto the continuum (|Dalgic et al 



2005). The HPQCD computa- 



1998b and uses 1-loop lattice 



2004|) . Therefore, they include 



all corrections to the heavy-light current through O (Aqcd / mb) > 0(a s ), 0(aa s ), O (a s / (amt,)) 
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and 0(a s AQCT>/ m b)- The HPQCD collaboration uses HMS5CPT for the chiral extrapolations 
of 4>b and 4>b s in a similar manner to Fermilab/MILC. They multiply the NLO expression by 
[1 +ca s a +c'a 4 ] in order to parameterize higher-order discretization effects. They also include 
an additional NNLO analytic term oc {m L i — m s ) in the extrapolation of the ratio 3>r /4>b. 



The HPQCD collaboration's ca 



20081) employs the HISQ action (IFollana et al 



culation of th e D an d D v -meson decay constants (|Follana et all 



2007|) (see Sec. III.EI) for all of the u, d, s, and 



c valence quarks. Because they are treating the charm quark as a light quark, the computation 
is similar to the determinations of f n and fa described in Sec. [Vlj, except for differences due to 
the fact that this is a mixed-action simulation with HISQ valence quarks and asqtad sea quarks. 
They use the medium-coarse, coarse, and fine MILC lattices, and include seven full QCD points 
in their analysis. They fix the c-quark mass so that the mass of the taste Goldstone r\ c meson 
agrees with experiment. Because the HISQ axial current is partially-conserved, it does not need 
to be renormalized. Therefore this method avoids the use of perturbation theory, whose truncation 
errors can be difficult to estimate. The HPQCD calculation does not use HMSXPT for the chiral 
extrapolations of fo and /o s , but simply applies continuum XPT, supplemented by Bayesian fit 
parameters. These parameters test for the expected discretization effects of the form asa 2 , a 4 , 
a? s a 2 , aga 2 log(m qimr j c ), an d 0? s a 2 m quar k from the asqtad action, and the effects of residual taste- 
violating interactions with HISQ valence quarks. 

All of the 2+1 flavor calculations of heavy-light meson decay constants rely upon power- 
counting in order to estimate the size of heavy-quark discretization errors. In the Fermilab method, 
heavy-quark discretization errors arise due to the short-distance mismatch of higher-dimension op- 
erators in the continuum and lattice theories. Th e sizes of these mis matches are estimated using 



HQET as a theory of cutoff effects, as described in 



Kronfeldl (12000) and 



Harada et al. 



(l2002bl) . This 



typically leads to errors of a few percent on the fine MILC lattices. In simulations with NRQCD b 
quarks, relativistic errors arise from higher-order corrections to the NRQCD action and the heavy- 
light current. Although these are not all discretization errors proportional to powers of the lattice 
spacing, many are proportional to inverse powers of the heavy-meson mass, and hence should 
be considered heavy-quark errors. The leading relativistic error comes from radiati ve corrections 
to the o B term in the action, and is estimated to be of 0{<x s Aqcq/ Mg) ~ 3% (|Gamiz et all 
2009|). The HISQ action is high ly-improved, a nd the leading heavy-quark errors are formally 



of o(a s (m c a) 2 ) and 0((m c a) 4 ) (|Follana et al. . 



2007|) . where a s ~ 0.3 and am c ~ 0.5 on the 
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fine MILC lattices. The HPQCD collaboration, however, removed errors of o(a s (m c d) ) in 
the HISQ action by accounting for radiative corrections in the coefficient of the Naik term, and 
also extended the traditional Symanzik analysis to remove all 0((m c a) 4 ) errors to leading order 
in the charm quark's velocity. Thus the leading charm quark discretization errors should be of 
((m c a) 4 (v / c) 2 ) ~ 0.5% or less for D mesons. 



C. Results for masses, decay constants, and CKM matrix elements 



Although the heavy-light meson decay constants, in combination with experimental measure- 
ments of leptonic branching fractions, can be used to extract CKM matrix elements via the relation 



T(H 



G 2 F \V ab \ 



flm 2 p M H I 1 
8ti JH 1 H \ M 2 



m 



(155) 



the matrix elements \V c d\, \V CS \, and \V u b\ can be obtaine d to better accuracy from other quanti 



2008). Therefore lattice 



ties such as neutrino scattering and semileptonic decays (lAmsler et all 
calculations of heavy-light meson decay constants provide good tests of lattice QCD methods, es- 
pecially the treatment of heavy quarks on the lattice. The comparison of lattice calculations with 
experimental measurements, however, relies upon the assumption that, because leptonic decays 
occur at tree-level in the standard model, they do not receive large corrections from new physics. 
This is generally true of most bey ond-the- standard m odel theories, but in a few mod els, such as 
those with leptoquarks, this is not necessarily the case (IDobrescu and KronfeldU2008|) . 

CKM unitarity implies that \V c d\ = \V US \ and \V CS \ = \V u d\ up to corrections of 0(|V U . S | 4 ). Because 
both \Vud\ and \V US \ are known to sub-percent accuracy, experimentalists use this relation to extract 
the D-meson decay consta nts from t he measured bran ching fractions. The latest determinations of 



fo (|Eisenstein et ail 120081) and fo s (|Alexandei , 



200% from the CLEO experiment are 
205.8 ± 8.9 MeV, f D + = 259.5 ±7.3 MeV . 



(156) 



These re sults use the determination of \V llc i\ = 0.97418(26 ) from superallowed 0+ — > + nuclear 



P-decay (|Towner and Hardy . 



2008) and of \V US \ = 0.2256 (lEisenstein et al. . 



2008 ). 12 The Fermi- 



lab Lattice and MILC collaborations' preliminary determination of the D-meson decay constants 



12 Although 



Eisenstein et al. 



(2008) attribute \V US \ = 0.2256 to FlaviaNet (lAntonelli . 



2007) 



Antonelli (2007) gives 



\V US \ = 0.2246(12) from decays plus lattice QCD, and |V„. S | = 0.2253(9) from K(2 and Kg^ decays plus lattice 
QCD. 
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FIG. 27 Comparison of lattice QCD and experimental results for fo and f Ds (left panel) and of lattice QCD 
results for fs and fg (right panel). 



are ( Bernard et a/.Ll2009b ) 



/d — 207(11) MeV, f Ds = 249(11) MeV 



(157) 



where the dominant errors come from tuning the charm quark mass and from heavy-quark dis- 
cretization effects, which are each ~ 3%. Both of these results are consistent with experiment. The 
HPQCD colla boration's determina tions of the D-meson decay constants using HISQ fermions are 



more precise (|Follana et all 



2007) 



f D = 207(4) MeV, f Ds = 241 (3) MeV 



(158) 



with total errors each below 2%. The largest contribution to the errors comes from the uncertainty 
in the scale r\, and is 1.4% (1%) for fo (/b s ). Although HPQCD's result for fo is consistent with 
experiment, their value for fo s is ~ 2.5-G below the CLEO measurement, where a is dominated 
by the experimental uncertainty. A comparison of lattice QCD and experimental results for the 
D-meson decay constants is shown in the left panel of Fig. |27] 

Many of the statistical and systematic uncertainties that enter the lattice calculations of fo and 
fr> s cancel in the ratio. Therefore the quantity /d/Id, allows for a more stringent comparison be- 



tween the results of Fermi 



find (Bernard etaL 



2009b) 



ab/MILC and HPQCD. The Fermilab Lattice and MILC collaborations 
fo/f Ds = 0.833(19), (159) 



while the HPQCD collaboration finds (|Follana et all 



2007): 



f D /f Ds = 0.859(8). 
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(160) 



The lattice results for the ratio disagree slightly, but only by ~ 1.6-G. The experimenta 



uncertain- 



tiesin fo and fr> s are largely independent, and therefore add in quadrature in the ratio ([Alexander , 
200% 

f D +/f D f= 0-793 ±0.040. (161) 
This increases the experimental errors and reduces the significance of the discrepancy with 
HPQCD. 

The HPQCD collaborati on also uses HISQ charm quarks to compute the D- and D v -meson 



masses (Follana et al. 



2007) 



M D = 1.868(7) GeV, M Ds = 1.962(6) GeV, 



(162) 



and their results agree with the experimental values Mq = 1.869 GeV and Mq s 



1.968 GeV (lAmsler et all 



2008). This lends credibility to their calculation of fo s , and suggests 



that both improved experimental measurements and lattice calculations are necessary to determine 
whether or not this discrepancy is new physics, a statistical fluctuation, or yet something else. 
Currently, Fermilab/MILC's determination of the -meson decay constant lies between the ex- 
perimental measurement and the calculation of HPQCD. Once the uncertainties in the calculation 
are reduced, which is expected to occur with the addition of statistics, finer lattice spacings, and a 
more sophisticated analysis, the Fermilab Lattice and MILC collaborations hope to shed light on 
this intriguing puzzle. 

5-meson leptonic decays are much more difficult to observe than D-meson decays because they 
are CKM suppressed (<=c |V J( b| 2 ). In addition, 5-decays to light leptons are suppressed by the factor 



mj in Eq. (11551) . and only decays to x's have been obser ved thus far, 



fraction T(B — > xv) is known only to ~ 30% accuracy (|Amsler et al. . 



urther more, the branching 



20081) . Thus there are no 



precise experimental determinations of the 5-meson decay constants, and the lattice calculations 
of fs and fs s should be considered predictions that have yet to be either confirmed or refuted by 
experiment. 

The F ermilab Lattice and M ILC collaborations preliminary determinations of fs, /b s , and the 



ratio are (Bernard et al. 



2009b) 



f B = 195(11) MeV, f Bs = 243(11) MeV, f B /f Bs = 0.803(28). 



(163) 



The largest errors in the individual decay constants are due to scale and light-quark mass uncer- 
tainties, light-quark discretization effects, and heavy-quark discretization effects, all of which are 
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~ 2 %. The HPQCD col laboration's determinations are consistent and have similar total uncertain- 



Their largest source of error is the ~ 4% uncertainty from 1-loop perturbative operator matching. 
A comparison of lattice QCD calculations of the 5-meson decay constants is shown in the right 
panel of Fig. [27] 

There are currently no calculations of the B- and B s -meson masses using the 2+1 flavor MILC 
lattices. This is, in part, because the staggered %PT expressions for heavy-light meson masses 
needed to extrapolate the numerical lattice data to the physical light-quark masses and the contin- 
uum are not known, and would require a nontrivial extension of the continuum expressions. 

MIL SEMILEPTONIC FORM FACTORS 

Lattice calculations of semileptonic form factors allow the extraction of many of the CKM 
matrix elements from experiment. The processes we consider for this purpose are dominated by 
tree-level weak decays of quarks at short distances, but are dressed by the strong interactions at 
longer distances, such that only mesons appear on the external legs. Given the nonperturbative 
form factor that parameterizes the strong interactions of the mesons, one can extract the CKM 
parameters that accompany the flavor-changing weak vertex. With enough processes one can over- 
constrain the four standard model parameters that appear in the CKM matrix, and thus test the 
standard model. 

A. D -> 7i£v and D -► K£v 

Semileptonic decays of D mesons, D — > Kiv and D — > n£v, allow determinations of the 
CKM matrix elements \V CS \ and \V c d\, respectively. Since these CKM matrix elements are well- 
determined within the standard model by unitarity, with results for other processes, the form factors 
can be obtained from experiment (assuming the standard model), and thus serve as a strong check 
of lattice calculations. Such calculations bolster confidence in similar calculations of B — > n£v, 
allowing a reliable determination of \V u b\, one of the more important constraints on new physics 
in the flavor sector. Precise calculations of semileptonic form factors for charm decays are also 




f B = 190(13) MeV, f Bs = 231(15) MeV, f B Jf Bs = 0.812(19). 



(164) 
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interesting in their own right, given the discrepancy between the HPQCD and experimental values 
for the D s leptonic decay. 

The necessary hadronic amplitude {P\V^\D) (P = K, it) is parameterized in terms of form factors 

by 

(P\V,\D) =f + (q 2 )(p D + pp-A), + fo(q 2 )A M , (165) 

where q = po~ Pp, \ = (mj) — m p)^/^ 2 > and V M = q^^Q. The differential decay rate dT/dq 2 is 
proportional to |V cx | 2 |/ + (^ 2 ) | 2 , with* = d, s. The CKM matrix element \V CX \ is determined using 
the experimental decay rate and the integral over q 2 of the lattice determination of \f + (q 2 )\. 

The matrix element (P|V^|D) is extracted from the three-point function, where the P meson is 
given a nonzero momentum p, 



C^ p ( W ,p) = ^y(O P (0)V,(y)ol(x)}, 
and Oo and Op are the interpolating operators for the initial and final meson st ates. The ca 



of this quantity by the Fermilab Lattice, MILC and HPQCD Collaborations (lAubin et al. 



(166) 



culation 



2005b) 



uses the Fermilab action [improved through 0(AQCD/ m c), with Aqcd in the HQET context] for 
the c quark and the asqtad action for the light valence quarks. The D meson and the heavy-light 
bilinears V M are constructed from a staggered light quar k an d a Wilson-type (Fe rmilab) heavy quark 



using the procedure described in lWingate et al 



(12003b and lBailevg? al. 



( 2009 ). In order to extract 



the transition amplitude (P|V^|D) from Eq. (11661) . we need the analogous two-point correlation 
function, 



Cffep) 



jyP- x <O w (0)Oi(x)) with M = D,P. 



(167) 



As in the case of decay constants, the renormalization factor matching the heavy-light currents on 
the lattice to the continuum is 

(168) 



Z$ q = p$ q Jz® Q z qq 

Vl_4 r Vi & V Va Va 



'Via 



J V 4 ^V 4 ' 



where the factors Zy® and Z qq are computed nonperturbatively, and th e remaining factor pffi 



2002b). 



(close to 1 by construction) is determined in one-loop perturbation theory (IHarada et all 

The quantities /y and f± are more natural quantities than / + and fo in the heavy-quark effective 
theory, and are defined as 



(P\V\D) = V^w[v"f\\(E)+plME)} : 



(169) 
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where v = PD/ m D, P_l = Pp — Ev and E = v ■ pp is the energy of the light meson. The chiral 
extrapolation and momentum extrapolation/interpolation are carried out in t erms of these param- 



Aubin et al. 



(2005b) 



eters, which are then converted into fo and / + . The chiral extrapolation in 
was performed at fixed E, wh ere /j| and f± were fit simultaneously to the parameterization of 
Becirevic and Kaidalovl (12000) (BK), 



Mq 2 



-, Hq 2 ) 



(170) 



(l-q 2 )(l-aq 2 Y ^ ' l-^/P' 
where q 2 = q 2 /m 2 ~ ) *, and F = f+(0), a and [3 are fit parameters. The BK form contains the pole 
in f+(q 2 ) at q 2 = m^. Even so, the BK parameterization builds into the calculation unnecessary 
model dependence. The more recent calculation of the similar semileptonic process B — > n£v does 
no t make use of this as sumption, as described in the next subsection. 



Aubin et al. 



(|2005bl) obtain for the form factors at q 2 = 

yJ-*(0) =0.64(3)(6), f^ K (0) = 0.73(3)(7). 



(171) 



where the first error is statistical, and the second is systematic. They also determine the shape 
dependence of the form factor as a function of q 2 . This is sho wn in Fig. [28l along with 



experimental data from the Belle Collaborat ion (lAbe et al 



tion. Taking the most recent CLEO results (|Ge et al 
f+^ K (0)\V cs \ = 0.744(7)(5) we obtain 



2005) that confirms their predic- 
20091) f+^ K (0)\V cd \ = 0.143(5) (2) and 



(172) 



\V cd \ =0.223(8)(3)(23), |V C „| = 1.019(10)(7)(106), 

where the first error is the (experimental) statistical error, the second is the (experimental) system- 
atic error, and the third is the total lattice error. If we use unitarity along with \ V m i \ and \ V US \ , then we 
can use the CLEO measurements to predict the form factors. We then obtain fP~* n (0) = 0.634(25) 
and ff^ K (0) = 0.764(9), in good agreement with the result in Eq. (11711) . Clearly, the lattice error 
still dominates the uncertainties. The largest errors in the lattice calculation are due to discretiza- 
tion errors and statistics. Improved calculations at finer lattice spacings and higher statistics are 
underway. 



B. £^7i£vand \V ub \ 

Comparison between theory and experiment for B — > n£v has been more troublesome than for 
other lattice calculations in CKM physics. Leptonic decays and BB mixing amplitudes are de- 
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FIG. 28 Comparison of the Fermilab/MILC/HPQCD lattice prediction for the normalized D — > K£v form 
factor (bands) with the subsequent Belle results (diamonds). The orange (dark gray) band is the 1 a error 
band from s tatistics, and the yellow (light gray) band is the 1 G band for all errors added in quadrature. 
Figure from lKronfeldl (120061 ) . 



scribed by a single parameter. The semileptonic decays B — > D^*'£v and K — > %£v can be described 
to high accuracy by a normalization and a slope. For B — > k£\, on the other hand, the form factors 
have a complicated q 2 dependence. Lattice data have covered only the low momentum, high q 2 
end of the pion momentum spectrum, and errors are highly g 2 -dependent and highly correlated 
between q 2 bins in both theory and experiment. 

It has long been understood that analyticity, unit arity, and cross i ng symmetry can be used 



to constrain the 



possib 



Boyd and Savagei 



1997 



e shapes of fo rm factors (|Bourrely et al. 



Lellouchi 



1981 



Boyd et all 



1995; 



1996|) . This has been used recently to simplify the compari- 
son of theory and experiment for B — > n£v. All form factors are analytic functions of q 2 except 
at physical poles and threshold branch points. In the case of the B — > %lv form factors, f(q 2 ) is 
analytic below the Btz production region except at the location of the B* pole. The fact that analytic 
functions can always be expressed as convergent power series allows the form factors to be written 
in a particularly useful manner. 
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Consider mapping the variable q 2 onto a new variable, z, in the following way: 



z(q 2 Jo) 



^\-q 2 /t+-^\-t /t + 
^l-q 2 /t + +^l-t /t+ 



(173) 



where t+ = (mg+mj) 2 , t- = (mg — m n ) 2 , and to is a free parameter. Although this mapping appears 
complicated, it actually has a simple interpretation in terms of q 2 ; this transformation maps q 2 > t + 
(the production region) onto \z\ = 1 and maps q 2 < t+ (which includes the semileptonic region) 
onto real z E [—1,1]. In the case of B — > n£v, the physical decay region is mapped into roughly 
—0.3 < z < 0.3. In terms of z, the form factors can be written in a simple form: 

1 



ftf) = n ( 2w 2 I S £ a k (t )z(q 2 M k - 



(174) 



Most of the q 2 dependence is contained in the first two, perturbatively calculable, factors. The 
Blaschke factor P(q 2 ) is a function that contains subthreshold poles and the outer function §(q 2 , to) 
is an arbitrary analytic function (outsid e the cut from t+ < q 2 < °°) which is ch osen to give the 



Bai ley et al. 



(2009), 



Arnesen et al 



(|2005|) . and references 



series coefficients a k a simple form. See 
therein for the explicit forms of these expressions. With the proper choice of §(q 2 ,to), analyticity 
and unitarity require the to satisfy 

N 

The fact that —0.3 < z < 0.3 means that according to analyticity and unitarity, only five or six 
terms are required to describe the form factors to 1% accuracy. (In B — > D^£v and K — > %£v, 
z is on the order of a few per cent in the physics decay region, which is why these decays can 
be accurately described by just two parameters.) Becher and Hill have argued that the heavy- 
quark ex pansion implies that the bound is actually much tighter than analyticity and unitarity alone 



demand (B echer and Hill, 



20061) . They argue that £f =0 a| should be of order (Aqcd/^) • This 



would lead to the expectation that only two or three terms will be sufficient to describe the form 
factors to 1% precision. 

Calculations have been performed by Fermilab Lattice and MILC collaborations using Fermi- 
lab b quarks, and by the HPQCD collaboration using NRQCD b quarks. Many of the details of the 
Fermilab/MILC calculations are the same as those for the Fermilab/MILC computation of heavy- 
light decay constants, described previously. For the semileptonic decays, only full QCD valence 
masses are used, as opposed to the partially-quenched masses used in leptonic decays. The calcu- 
lations use the a ~ 0.12 and 0.09 fm gauge field ensembles. The HMSXPT continuum and chiral 
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FIG. 29 Results for the normalized B — ► n£v form factor P 



lations (circles) and BABAR (stars), from 



Bailey et al. 



(j) + / + from the Fermilab/MILC lattice calcu- 



((2009J). The solid (red) line is the results of a fully 



correlated simultaneous fit. Requiring that lattice and experiment have the same normalization yields |V i( fe|- 

extrapolations are done with the full NLO expressions plus additional NNLO analytic terms. These 
formulae allow the simultaneous interpolation in pion energy along with the continuum and chiral 
extrapolations, thus reducing the total systematic uncertainty. 

Figure [29] shows the result of a fully correlate d simultaneous z-fit to the Fermilab/MILC lattice 



data and the BABAR 12-bin experimental results (|Aubert et a/.U2007|) . with \ V u b \ being a parameter 
in the fit. The resulting z-fit parameters are a = 0.0218 ± 0.0021, a\ = -0.0301 ± 0.0063, a 2 
-0.059 ± 0.032, a 3 = 0.079 ± 0.068, and 



\Vub\ = (3.38 ±0.36) x 10 



-3 



(176) 



(Bai ley et al 



2009 ). The coefficients of z" are indeed of order (AQco/ m &) 3//2 a s argued by 



Becher and Hilll (120061) . Because the ~ 11% uncertainty comes from a simultaneous fit of the 
lattice and experimental data, it contains both the experimental and theoretical errors in a way that 
is not simple to disentangle. If we make the assumption that the error in \V u t\ is dominated by the 
most precisely determined lattice point, we can estimate that the contributions are roughly equally 
divided as ~ 6% lattice statistical and chiral extrapolation (combined), ~ 6% lattice systematic, and 
~ 6% experimental. The largest lattice systematic uncertainties are heavy-quark discretization, the 
perturbative correction, and the uncertainty in gB*Bn, all of which are about 3%. Our determination 
is ~ 1 — 2o lower than most inc l usive determinations of | V„b\, where the values tend to range from 



4.0-4.5 x 10 3 ( Pi Lodovico . 



20081) . Our determination is, however, in good agreement with 
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the pr eferred values from the CKMfitter Collaboration (\V llb \ = (3. 44 + gffi x 
20081) ) and the UTfit Collaboration (\V ub \ = (3.48 ±0.16) x 10~ 3 Jsilvestrini . 



3 ( Charles et all 



2008)). 



Many of the details of the HPQCD calculation of B — > n£v are the same as described for heavy- 
light decay constants in the previous section. They use NRQCD b quarks and asqtad light quarks. 
On the coarse, a « 0.12 fm ensembles, they perform the calculation on four unquenched ensem- 
bles plus an additional two partially quenched light quark masses on one ensemble. They also use 
full QCD data on two fine, a ~ 0.09 fm ensembles in order to constrain the size of discretization 
effects. They use HMSXPT to perform the chiral extrapolations separately for various fiducial 
values of E % after interpolating in E n . They also show that they obtain consistent results with sim- 
pler chiral extrapolation methods. They perform fits to their data using the z-fit method described 
abov e, as well as several oth e r func tional forms including t he Becirevic-Kaidalov p arameteriza- 



tion (Becirevic and Kaidalov, 



20001) and Ball-Zwicky form (IB all and Zwicky , 



2005). Note that 



they do not use a combined fit of experimental and lattice data using the z-fit method to extract 
\Vub\. Rather, they use the various parameterizations to integrate the form factor f+(q 2 ) over q 2 , 
and they show that they obtain consistent results wi th all methods. App lying their results to 2008 



data from Heavy Flavor Averaging Group (HFAG) (|Di Lodovicol |2008|) yields 



|V^| = (3.40±0.20+°;g)xlO 



(177) 



(IDalgic et all 120061) . where the first error is experimental and the second is from the lattice calcu- 



lation. Figure [30] shows the comparison between an average of the Fermilab/MILC and HPQCD 



results for \V u f,\ and two inc 



inputs (IGambino et al. 



2007 



usive determination s of the same quantity using different theoretical 



Lange et all 



2005). 



C. B -> Dlv and B -> D*N 

The CKM parameter \V c b\ is important because it normalizes the unitarity triangle character- 
izing CP-violation in the standard model, and must be determined precisely in order to constrain 
new physics in the flavor sector. The standard model prediction for kaon mixing contains \V c b\ 
to the fourth power, for example. It is possibl e to obtain JV^ I f r om both inc l usive and exclusive 



semi 



eptonic B decays. The inclusive decays (IBigi et al. . 



1992a . 



1997 



1993, 



1992b; 



Chay et all 



1990|) make use of the heavy-quark expansion and perturbation theory, while the exclusive decays 



require the lattice calculation of the relevant form-factors. Each of the exclusive channels B — > D£v 
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ub 

FIG. 30 Values of \V u b\ obtained from averaging the exclusive determinations compared with inclusive 
determinations using different theoretical inputs. 

and B — > D*£v allows a lattice extraction of \V c b\, and thus they provide a useful cross-check, both 
of each other, and of the inclusive determination. We have so far considered the calculations of the 



necessary form fa ctors only at zero-recoi 



theoretical errors (Hashi moto et all 



2002). 



, as this leads to considerable simplification and reduced 



The differential rate for the decay B — > D£\ is 
dT(B -> D£v) Gj 



dw 



48tc 3 



with 



g(w) 



m 3 D (m B + m D ) 2 (w 2 -\) 3 / 2 \V cb \ 2 \g(w)\ : 



m B -m D 
h+(w) ; h-{w), 



(178) 



(179) 



m B + m D 

where Gf is Fermi's constant, h + (w) and h-(w) are form factors, and w = V ■ v is the velocity 
transfer from the initial state to the final state. The differential rate for the semileptonic decay 
B -> D*£v e is 



dT(B^D*£v) Gl , . n2 /-= . l2 . ,. . . 

rfitjy* (m B - m D *) 2 Vw 2 - 1 \V cb \ 2 %(w) 1 7 (w) 



(180) 



dw 47l 3 

where %(w)|5 (w)\ 2 contains a combination of four form factors that must be calculated nonper- 
turbatively. At zero recoil (w = 1) we have %(1) = 1, and f (1) reduces to a single form factor, 

Mi). 



We compute the form factor h + at zero-recoil using the double ratio (|Hashimoto et all 



1999) 



(D\cY 4 b\B)(B\by 4 c\D) _ 2 



(181) 



(D\cy 4 c\D)(B\by 4 b\B) 

This double ratio has the advantage that the statistical errors and many of the systematic er- 
rors cancel. The discretization errors are suppressed by inverse powers of heavy-quark mass as 
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a^(Agci>/2mg) 2 and (Agcz?/2mg) 3 (|KronfeldL l200Qh . and much of the c urrent renormalizatio n 
cancels, leaving only a small correction that can be computed perturbatively (|Harada et a/.Ll2002ah 
The extra suppression of discretization errors by a factor of A/ lmg occurs a t zero-recoil for heavy- 



1990). 



to-heavy transitions, and is a consequence of Luke's Theorem (|Lukd. 

In order to obtain /z_, it is necessary to consider nonzero recoil momenta. In this case, Luke's 
theorem does not apply, and the HQET power counting leads to larger heavy-quark discretization 
errors. However, this is mitigated by the small con tribution of /?_ to the br anching fraction. The 



form factor h- is determined from the double ratio ([Hashimoto et all 

(D\cjjb\B) (D\cy 4 c\D) 



1999) 



(D\cy 4 b\B)(D\cjjb\D) 



1 



h-{w) 
h+(w) 



2h+ (w) 



(182) 



which is extrapolated to the zero-recoil point w = 1. Combining the determinatio ns of h + (l) and 



/?_(!), we obtain the preliminary result (?(1) = 1.082(18) (16) (|Okamoto 



2006J), where the first 



error is statistical and the second is the s um of all sys tematic errors in quadrature, and where 



982) . Combining this with the latest average 



we have included a 0.7% QED correction (ISirlinl. 

from the (HFAG), Q {l)\V cb \ = (42.3±1.5) x 10~ 3 (iDi Lodovicol l2008l). we obtain the preliminary 
result 



\V cb \ = (39.1±1.4±0.9) x 10 



-3 



(183) 



where the first error is experimental, and the second is theoretical. 

The form factor at zero-recoil needed for B — > D*£v is computed using the double ratio 



( Bernard et al. 



2009a) 



{D*\cyjy 5 b\B)(B\byjy 5 c\D*) 



imi; 



(184) 



(D*\cj 4 c\D*)(B\by 4 b\B) 

where again, the discretization errors are suppressed by inverse powers of heavy-quark mass as 
a 5 (A,2C£)/2mg) 2 and (Agcz)/2mg) 3 , and much of the cur rent renormalization cancels, leaving 



only a small correction that can be computed perturbatively (IHarada et al. 



to physical light quark masses using the appr opriate r HMSX PT (|Laiho and Van de Waten . 



({Bernard et al. . 



Including a QED correction of 0.7% (ISirlinl Il982l) . we obtain j(l) = 0.927(13)(20) 



2002 aj). We extrapo 



2006). 



ate 



2009al) . where the first error is statistical and the second is the sum of system- 



atic errors in quadrature. Taking the latest HFAG av erage of the experimental determination 
7 (l)|Vd,| = (35.49 ±0.48) x 10~ 3 JpiLodovico . 



20081) . we obtain 
\V cb \ = (38.3 ±0.5 ±1.0) x 10~ 3 , 



(185) 
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FIG. 31 Values of \V c b\ from the exclusive decays B — > Dlv, B — > D*£v, and the inclusive determination. 



The experimental average includes all available measurements of J (1) |V C ^|, but we point out that 
the global fit is not very consistent [% 2 /dof = 39/21 (CL=0.01%)]. The Particle Data Group 
handles this inconsistency by inflating the experimental error by 50% (lAmsler et all 120081) . The 



dominant lattice errors are discretization errors and statistics, and work is in progress to reduce 
these. Note tha t there is some tensio n between this and the inclusive determination of \V c b\ = 
41.6(6) x 10~ 3 dBarberio et a/1120071) , as can be seen in Fig. 



IX. OTHER COMPUTATIONS USING MILC LATTICES 



In this section, we describe a variety of additional results based on the MILC ensembles. Over 
eighty-five physicists outside the MILC collaboration have used the MILC configurations in their 
research. This includes colleagues at nearly forty institutions throughout the world. Their research 
covers a very broad range of topics including determinations of the strong coupling constant, the 
quark masses, the quarkonium spectrum and decay widths, the mass spectrum of mesons with a 
heavy quark and a light antiquark, the masses of baryons with one or more heavy quarks, as well 
as studies of the weak decays of mesons containing heavy quarks, the mixing of neutral K and B 
mesons with their antiparticles, the quark and gluon structure of hadrons, the scattering lengths of 
pions, kaons and nucleons, the hadronic contributions to the muon anomalous magnetic moment, 
and meson spectral functions. 
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A. Determination of the strong coupling constant and the charm quark mass 



1. The strong coupling constant from small Wilson loops 



The HPQCD collaboration used MILC 



stant a s ( Davies et all 



2008 



Mason et all 



attice ensembles to compute the strong coupling con- 



2005b . They compute nonperturbatively (i.e., numeri- 



cally on the MILC lattices) a variety of short-distance quantities Y, each of which has a perturbative 
expansion of the form 

oo 

Y = £c n c$(d/a) , (186) 

n=l 

where c„ and d are dimensionless q-indepe ndent constants, and av / fJ/a) is the running QCD 
coupling constant in the so-called V scheme (|Lepage and Mackenzie! . 1 1993b for rif = 3 flavors of 
light quarks. 

The coupling ay(d/a) is determined by matching the perturbative expansion, Eq. (11861) . to the 
nonperturbative value for Y. Perturbatively converting from the V to the MS scheme and running 
up to the Z boson mass, switching to rif = 4 and then 5 at the c and b quark masses, gives a 
determination of the strong coupling constant CX^(Mz,n/ = 5). 

The HPQCD collaboration considered 22 short distance qu antities Y, consistin g of the loga- 



rithms of small Wilson loops and ratios of small Wilson lo ops (IDavies et al. 



in Eq. (11861) are determined perturbatively by the method of [Lepage and Mackenzie! (| 1993b . c n for 



2008). The scales d 



n = 1, 2 and 3 were computed in lattice perturbation theory (|MasonL 12004b . and higher orders, up 
to n = 10 were included in a constrained fitting procedure. In practice, dy(d/a) for all the different 
scales d/a used was run to a common scale of 7.5 GeV, and oto = ay (7.5 GeV) was used as a free 
fitting parameter in the constrained fits for each of the observables. 

Corrections to the perturbative form, Eq. (11861) . from condensates appearing in an operator 
product expansion (OPE) for short-distance objects, were included in the constrained fitting pro- 
cedure. Other systematic errors such as finite lattice spacing effects and scale-setting uncertainties 
were considered. As their final result, the HPQCD collaboration quotes 

av/(7.5GeV,n/ = 3) =0.2120(28) and a m (M z ,n f = 5) = 0.1183(8) . (187) 

The l attice determina t ion of a^? (M z ) is compared to other determinations in Fig. 



In 



Maltman et al. 



(|2008b a reanalysis of three of the short distance quantities used by the 
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FIG. 32 Summary of determinations of the strong coupling constant a s (Mz) from lAmsler et al\ (I2008r) . The 
lattice QCD determination is the most precise one. 



HPQCD collaboration was performed with the result 



a m (M z ,n f = 5) =0.1192(11) 



(188) 



2007). The two 



in good agreement with other next-next-to-leading-order determinations (Bethke, 
analyses differ in the way the perturbative running and matching was done, the value of the 
gluon condensate used in the OPE subtraction, the way the scale setting for each lattice ensem- 
ble is treated and a sl ight difference of the value used for the scale setting. For more details see 



Maltman et al. 



(2008). 
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2. The charm quark mass and the strong coupling constant from current-current correlators 



A new approach to extract a s and to determine the charm quark mass m c was used in 



Allison et al. 



(|2008l) . It consists of comparing moments of charmonium current-current correlators 



computed nonperturbatively on the lattice with high-order continuum QCD perturbation theory. 
Vector current-current correlators have previously been used to obtain so me of the most precise 



determinations of m r fro m the experimental e + e — »■ hadrons cross section ( 



2001 



Kiihn et al. . 



Kiihn and Steinhausei . 



2007). On the lattice, many types of correlators are available that are not acces- 
sible to experiment. In particular, the pseudoscalar current-current correlator can be computed to 
very high statistical accuracy, and the presence of a partially-conserved axial vector current makes 
current renormalization unnecessary. 

Consider the r\ c current-current correlator 

G{t) =a 6 £(am , c ) 2 (0|j 5 (x,0j5(0,0)|0) , (189) 

X 

with moments 

T/2 

G n = £ (t/a) n G(t). (190) 

t=-T/2 

In the continuum limit, these moments can be computed perturbatively as 



gn («Ms(//),/ //m c ) 
(am c (/j)) r 



R, = G A /Gf and R n = ^ (gJG^) ^ for n > 4 . (192) 



where g n is known to O (oc^) for n = 4, 6 and 8. The approach to the continuum limit is improved 
by dividing by the tree-level results, and tuning errors in m c and errors in the scale setting are 
ameliorated by multiplying with factors of the lattice T| c mass 

2am 0c 

The ratios R n are extrapolated to the continuum limit using constrained fits. Comparing with 
continuum perturbative ratios r\ = g^/gf^ and r n = (g n /g^) 1 ^' 1 for n > 4, allows a^g to be 
extracted from R4 and ratios R n /R n+ 2 given the charm quark mass, and the charm quark mass can 
be obtained from the R n with n > 4, given the value of the strong coupling constant, 

cW 2 R n (a = 0) • ( 
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Allison et al.\ (120081) used eight MILC lattice ensembles with fo ur different lattice spacing s. The 



charm correlators were computed using HISQ staggered quarks (|Follana et al. 
obtained for m r 



2008, 



20071) . They 



m c (3GeV,n/ = 4) =0.986 (10) GeV , or m c (m c ,n f = 4) = 1.268(9) GeV. 



(194) 



his is in good a greement, and about twice as precise as the best previous determination 



(Kuhn et al. 



20071) . They obtain for a s 



a m (M z ,n f = 5) =0.1174(12) 



(195) 



in good agreement wit h the lattice determination described earlier and with other NNLO determi- 



nations (Bethke, 



2007) 



B. Onia and other heavy mesons 

Heavy quarkonia were important in the early days of QCD because potential models could be 
used to approximately understand their dynamics before first-principles calculations were possi- 
ble. The approximate validity of potential models helps in the selection of operators needed in the 
improvement program for quarkonia. The several methods for formulating heavy quarks on the 
lattice have various advantages and disadvantages for quarkonia. NRQCD employs the operators 
of the nonrelativistic, heavy-quark expansion. The operator expansion converges poorly for char- 
monium, and fails when Aqcd/^ is not small. The Fermilab interpretation of Wilson fermions 
interpolates between a nonrelativistic type of action at ma 1 and the usual Wilson-type action 
at ma <C 1. It can be used for all ma, but has a more cumbersome set of operators, and has been 
less highly improved than other heavy-quark actions. The HISQ action is a light quark action that 
fails when ma 1 , but has been improved at tree level to high orders in ma and works well for ma 
close to 1 . 



1. Bottomonium with NRQCD heavy quarks 



The HPQCD and UKQCD collaborations have studied bottomon ium spectroscopy on several 



MILC ensembles with lattice spacings a ~ 0.18, 0.12 and 0.09 fm (|Gray et al. 



2005 ). Even on 



the finest of these ensembles, am^ ~ 2. The authors have used lattice NRQCD to formulate the b 
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quarks in the regime am > 1 (|Davies et al. 



1994; 



Lepage et al 



1992; 



Tha cker and Lep age. 



1991) 



The form of the action of NRQCD was shown in Eq. (11441) . The b quark is nonrelativistic inside 
the bottomonium bound states, with velocity v? ~ 0. 1 . NRQCD, as an effective field theory, can 
be matched order by order to full QCD in an expansion in v 2 and a s . The action currently in use 
includes corrections of (v 2 ) beyond leading order. Discretization errors have also been corrected 
to the same order in v 2 . 

The spin-averaged T mass splittings are expected to be quite insensitive to many lattice un- 
certainties, such as light sea quark masses and normalization of cor rection operators. They are, 



Gray et al. 



(|2005|) compute 



therefore, expected to be calculable to high accuracy on the lattice, 
spin-averaged mass splittings, IP — 15 (i.e., 1 — l 3 5i), 25 — 15 (i.e., 2 3 5i — l 3 5i), IP — IS, and 
35 — 15 in lattice units, and then use the experimental splittings to determine the lattice scale, as 
described in Sec. IIV.Bl Figure [33] shows the results, where the lattice spacing has been set by the 
25 — 15 splitting, and has been set from My. The left-hand figure compares the results in GeV 
at two lattice spacings, for quenched and unquenched calculations. The right-hand figures show 
the splittings calculated on the lattice divided by experiment, in the quenched approximation (left 
narrow figure) and unquenched (right narrow figure). Clear disagreements with experiment in the 
quenched approximation are removed in the unquenched calculations. 



As for the T(15) hyperfine splitting, 



Gray et al. 



(|2005|) quote AM = 61(14) MeV, corresponding 



to r\tM = 0.099(22), following an extrapolation to the physical point. This result is cons istent 



with the recent observation of the r\i, by the BABAR collaboration (|Aubert et al. . 
found a splitting of 71(4) MeV from the T(15). 



2008 



2009) who 



2. Onia with Fermilab quarks 

The Fermilab and MILC collaborations have computed charmonium and bottomonium masses 



( Gottlieb et al. . 


2006a 


4 


di Pierro et al. , 


2004) 



Fermilab quar ks (lEl-Khad 
In Fig.l34l (lDeTaref al. 



sl et al 



2OO9I all the resulting masses for charmonium and bottomonium are 



20041) . For the heavy charm and bottom quarks the y use 



19971) . An updated study is underway (|DeTar et all 



2009). 



shown as splittings from the spin-averaged 15 state. Plotted are the chirally-extrapolated values for 
each lattice spacing. They are compared with the experimental values given by solid lines, where 
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FIG. 33 Left: the T spectrum of spin-averaged radial and orbital levels in GeV. Closed and open symbols are 
from coarse and fine lattices respectively. Squares and triangles denote unquenched and quenched results, 
respectively. Lines represent experiment. Right: spin-averaged mass differences from the same data divided 
by e xperiment, in the q uenched approximation (left narrow figure) and unquenched (right narrow figure), 



from 



Gray et al. 



(2005). 
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FIG. 34 Summary of the charmonium (left) and bottomonium (right) spectra. The fine ensemble results 
are in blue fancy squares, the coarse in green circles, the medium coarse are in orange diamonds and the 
extracoarse results are in red squares. Included in the error budget is an estimated systematic uncertainty 
from setting the heavy quark masses. 

the experimental results are known. In the cases where they are not known and are estimated from 
potential models, they are shown as dashed lines. The charmonium spectrum shows good agree- 
ment with experiment for the ground states, except for the % c o, which may be slightly heavier than 
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the experimentally measured value. The excited S- wave states are also heavier than their respective 
experimental results, but one has to bear in mind that these states are difficult to determine without 
careful consideration of finite-volume effects since they are close to the DD threshold. The bot- 
tomonium summary panel shows the general tendency of the result to approach the experimental 
values as the lattice spacing decreases. 

Charm annihilation processes give a po ssible additional corr e ction to the charmonium hyperfine 



splitting. 



DeTar and Levkoval ( 20071) and 



Levkova and DeTarl (|2009h have started to study these 



quark-line disconnected diagrams using MILC ensembles wi th lattice spacings a ~ 06 and 0.09 
fm. They use stochastic estimators with unbiased subtraction (|Mathur and DongL 120031) to compute 
the disconnected contribution to the r| c propagator. They find that annihilation processes increase 
the r\ c mass a small amount (by 5.5(8) MeV for a fine lattice and 3.4(3) MeV for superfine), thereby 



decreasing slightly the predicted hyperfine splitting (|Levkova and DeTai , 



2009 



3. Charmonium with highly improved staggered quarks 



The HPQCD and UKQCD collaborations have studied charmonium spectroscopy on MILC en- 
sembles using the HISQ action for the valence quarks. They use MILC ensembles with lattice 
spacing a ~ 0.12 and 0.09 fm, where am c = 0.66 and 0.43, respectively, to demonstrate the advan- 
tages of the HISQ action, and compute the charmonium spectrum, using the T| c . mass to tune the 
input value for am c . They have corrected discretization errors in am up to order (am) 4 , and shown 
that this produces a speed of light that is independent of p a nd equal to 1, within errors, in the 



Follana et al. 



(|2007l) . In particular, 



equation E 2 = p 2 c 2 + m 2 c 4 . The results are shown in Fig. 7 of ] 
they find for the hyperfine mass splitting Mj/^ — M^ c = 109(5) MeV. This result is the closest to 
the physical value of 1 17(1) MeV that has yet been achieved. 



4. The B r meson 



The HPQCD, Fermila b Lattice and UKQC D collaborations used MILC ensembles to predict 



the mass of the B ( meson (Allison et al. 



2005|) before it was accurately measured. They used two 



different fermion actions for the heavy bottom and charm valence quarks, ch oosing the more op 



timal action in each case. For the bottom quark, they used lattice NRQCD (|Davies et all 



1994; 
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Lepage et all 119921 : iThacker and Lepagd . 119911) . because it has a better treatment of the v 4 inter- 
actions, where v is the velocity of the heavy quark. For th e charm quark, they used the relativistic 



Fermilab action (El-Kh adra et al. 



1997 : 



Kronfeldl 



200oL which treats higher order effects in v 2 



be tter. This is appropri ate, since the velocity of the c quark in B c is not particularly small, v c ~ 0.5. 



Allison et al. 



(|2005|) calculated mass splittings, for which many of the systematic errors cancel, 
namely 

A ¥ t = m Bc - (7% + m T ) /2 , A DsBs = m Bc - (m Ds + m Bs ) , (196) 

where m ¥ = (m Tlc + 3m y / v )/4, m Ds = (m Ds + 3m D *)/4, and m Bs = (m Bs +3m B *)/4 are spin- 
averaged masses. They found no visible lattice-spacing dependence using ensembles with a fa 
0. 18, 0. 12 and 0.09 fm. Extrapolating the a ~ 0. 12 fm results linearly in the light sea quark mass 
they obtain 

A ¥Y = 39.8±3.8±11.2+J 8 MeV, A DsBs = -[1238 ± 30 ± lljj[ 7 ]MeV • (197) 
The errors are from statistics, tuning of the heavy-quark masse s, and heavy-quark d iscretization 



effects. Since the statistical error on the first splitting is smaller, 
predict the B c mass as 

m Bc = 6304 ± 4 ± 1 1 +J 8 MeV . 



Allison et al. 



(|2005h used that to 



(198) 



Shortly after the lat tice calculation was pub lished, the CDF collaboration announced their precise 



mass measurement (lAbulencia et al. 



2006) 



m Br = 6287 ±5 MeV 



(199) 



in good agreement with the lattice prediction, i.e., slightly more than 1-a away. 



C. Heavy baryons 



Baryons containing a hea yy quark comprise a r ich set of states. For example, there are currently 
17 known charmed baryons (|Amsler et ail |2008|) . However, for bottom baryons, there are only a 
few known states. Thus, it is possible both to verify calculations by comparison with known masses 
and to make predictions for as yet undiscovered states. 

Many of the heavy baryons contain one or more u or d qua rks, thus requiring a chiral ex 



trapolation. Although some early work on MILC configurations (IGottlieb and Tamhankar 



2003; 
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Tam hankarl . 120021) used clover quarks for u, d and s, this limited how clos ely one could approach 



the c 



liral li mit, and recent work has used staggered light quarks instead (|Na and Gottlieb 



2009, 



2006 



20071) . The heavy quark is dealt with as in Sec. IVII.AI 



Th e pioneering lattice work on heavy baryons by the UKQCD collaboration (|Bowler et all 
19961) considered two operators 5 = £abc(tf\ CJs^^h and °n = ZabciY^Cy^)^, where 



e a t, c is the Levi-Civita tensor, \f/i and ^2 are light valence quark fields for up, down, or strange 
quarks, is the heavy valence quark field for the charm or the bottom quark, C is the charge 
conjugation matrix, and a, b, and c are color indices. The former operator can be used to study the 
spin- 1/2 baryons A/ ; and S/,. The latter can be used, in principle, for both spin- 1/2 and spin-3/2 
baryons. However, with the current formalism, for operators with two staggered quarks, there are 



cance llati ons in the spin - 3/2 se ctor and can only be used for spin- 1/2 baryons (|Na and Gottlieb 



20071) . In 



Gottlieb et 



al. 



in much the way that 



(120081) the taste properties of staggered di-quark operators are considered 



Baileyl (I2007|) studied staggered baryon operators. However, this method has 



not yet been applied in calculations. For states with two heavy quarks, both spin- 1/2 and spin-3/2 
states have been studied. 

Another issue when dealing with states containing heavy quarks is the distinction between the 
rest and kinetic masses (see Sec. I VIII) . Calculation of kinetic masses requires looking at states with 
nonzero momentum and fitting a dispersion relation. This has not yet been done for the heavy 
baryons, which means that we are restricted to reporting mass splitt ings. 



2009). With 



So far, ensembles with three lattice spacings have been studied (|Na and GottliebL 
a « 0. 15 fm, three ensembles with mi/m s = 0.2, 0.4 and 0.6 were used. With a « 0. 12 fm, mi/m s = 
0.007, 0.01 and 0.02, and with a ~ 0.09 fm, only m//m ? = 0.2 and 0.4 were studied. Seven to nine 
light quark masses are used to allow for chiral extrapolation. The charm and bottom quark masses 
are as in the meson work. Since mass splittings are desired, ratios of hadron propagators are fit 
in preference to fitting each hadron and subtracting the masses. For baryons with a heavy quark, 
rSZPT has not been worked out yet, so the chiral extrapolation is based on a polynomial in the 
valence and sea masses, 



^quad = c + c\mi + c 2 m, + c 3 m s + c 4 m sea 



(200) 



where cq to C4 are the fitting parameters, m/ is the light valence quark mass, m s is the strange 
valence quark mass, and m sea is the light sea quark mass. These fits are denoted "quad" in the 
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FIG. 35 Indepen dent mass differences o f J p = A + singly charmed baryons (a), and singly bottom baryons 



(b). Figures from 



Na and Gottlieb! (12009n . 



figures. Alternative chiral extrapolations use only the full QCD points, i.e., those in which the 
valence and sea light quark masses are equal. These are denoted "full" in the figures. 

For the singly-charmed baryons in Fig. l35fa). three of the four differences are in good agree- 
ment with the experimental results. The result that is not in good agreement is one that involves 
one hadron from O5 and one from O^. The other differences come from particles that are both 
determined using the same operator. This behavior is a mystery. 

In Fig. [35jb), we consider the singly-bottom baryons and find good agreement for the one 
observed difference for Eh — A/ 7 . Also shown is the comparison with a recent lattice calculation of 



Lewis and Woloshvn 



(120091) . The large value for the spli tting is again noticeable. 



In Fig. [361 we compare with the results of 



Lewis et al. 



(|200lh and 



Lewis and Woloshynl (12009) 



for both spin- 1/2 and spin-3/2 baryons. The earlier calculation of charmed baryons used quenched 
anisotropic lattices generated with an improved gauge action. The more recent calculation of bot- 
tom baryons uses configurations containing the effects of dynamical quarks. In order to compare 
the two calculations, and because kinetic masses are not available in the calculation on MILC con- 
figurations, a constant was added to the static masses that depends on lattice spacing and whether 
the state contains charm or bottom quarks, but not upon spin or light quark content. 

There are a number of ways to improve upon the current work including increasing statistics, 
extending the calculations to the finer ensembles, studying the kinetic masses and studying new 
operators that will allow us to explore the properties of the spin-3/2 baryons. It is also possible to 
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FIG. 36 The mass spectrum of dou bly charmed and bottom baryons. The error bars are statistical only. 



Figures from 



Na and Gottlieb! (I2009h . 



use HISQ quarks for all of u, d, s and c quarks to explore the charm sector using only staggered 
operators. 



D. K° - K° mixing: B K 

Experimental measurements of the size of indirect CP- violation in the neutral kaon sys tem eft- 
can b e combined with theoretical input to constrain the apex of the CKM unitarity triangle (IBuras 



2008), the 



1998). Because Zk has been measured to better than a percent accuracy (|Amsler et all 
dominant sources of error in this procedure are the theoretical uncertainties in the CKM matrix 
element \V c b\, which enters the constraint as the fourth power, and in the lattice determination of 
the nonperturbative constant Bk- 

The kaon bag-parameter Bk en codes the hadronic contribution to K° — Z° mix- 



ing (|Buchalla et ail 



1996 



Buras 



1998) 



'^°|£as=2(a<)I* > 



Bk(m) 



f(^|syoY 5 J|0>(0|JYoY5^ > 



(201) 



where Qas=2 is the effective weak four-fermion operator 

Qas=2 (x) = [fyd] v-a (x) [sJnd] v-a (x) (202) 

and /j is a renormalization scale. The dependence on /j cancels that of a Wilson coefficient C(/j) 
that multiplies Bk(h) in physical observables such as the mass difference between Ks and K^. 
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The denominator in Eq. (12011) is the value of the matrix element with vacuum saturation of the 
intermediate state. Often quoted is the value of the renormalization group invariant form of Bk, 
Bk, defined by 

B K = C( M )B K (n) . (203) 



Gamiz et al 



(|2006l) carried out a calculation of Bk using two MILC ensembles with lattice 
spacing a ~ 0. 12 fm. They employed asqtad valence quarks with valence kaons made of degenerate 
quarks of mass m s /2. Using one-loop matching with the coupling taken as ocy (1/a) they find the 



following value for Bk in the naive dimensional regularization scheme: 



B 



MS- NDR 
K 



(2GeV) =0.618(18)(19)(30)(130) 



(204) 



where the errors are from statistics, the chiral extrapolation (I Van de Water and Sharpel . l2006|) . dis- 
cretization errors, and the perturbative conversion to the MS — NDR scheme. The value Eq. (12041) 
corresponds to Bk = 0-83 ±0.18. The error is dominated by the uncertainty from O(a^) correc- 
tions to the perturbative lattice-to-continuum matching. 

Because the matching coefficients are known only to one loop, the result in Eq. (12041) is not 
competitive with the published domain-wall fermion calculation by the RBC and UKQCD Col- 
laborations, in which the operator renormali zatio is done nonperturbatively using the method of 
Rome-Southampton (IMartinelli et all 1 1995b and mixing is suppressed due to the approximate 



chiral symmetry. 



(Allton et al 



hey obtain, using a single, comparable lattice spacing, Bk = 0.720 ±0.019 
2008J), where the dominant uncertainty is due to discretization errors, and is esti- 
mated to be ~ 4% from the scaling behavior of quenched data. 

Recently Aubin, Laiho, and Van de Water obtained the first unquenched determi nation of B k 



at two lattice spacings using domain-wall valence quarks on the MILC ensembles (|Aubin et al 



2009bl) . Because dynamical domain-wall lattice simulations are computationally expensive, this 
mixed-action approach is an affordable compromise that takes advantage of the best properties of 
both fermion formulations. Since the MILC ensembles are available at several lattice spacings 
with light pion masses and large physical volumes, this allows for good control of the chiral ex- 
trapolation in the sea sector and the continuum extrapolation. Domain-wall fermions do not carry 
taste quantum numbers, so there is no mixing with operators of other tastes. Furthermore, the ap- 
proximate chiral symmetry of domain-wall fermions suppresses the mixing with wrong-chirality 
operators and allows the use of nonperturbative renormalization in the same manner as in the purely 
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domain- wall case. Finally, the e xpression for Bk in m ixed action %PT contains only two more pa 



rameters than in continuum %PT (Aubin et al 



2007b|) . both of which are known and are, therefore, 



not free parameters in the chiral and continuum extrapolation. 



Aubin et al. 



( 2009bl) obtain 



B 



MS- NDR 
K 



(2GeV) =0.527 (6) (20), 



(205) 



where the first error is statistical an d the second is systematic. With data on the coarse and fine 



MILC lattices, 



Aubin et al. 



(|2009br) find that the discretization errors in Bk are small. The largest 
error in Bk is ~ 3% and is from the renormalization factor Zg K , which is computed nonperturba- 
tively in the RI/MOM scheme, but must be converted to the MS-scheme using 1-loop continuum 
pe rturbation t heory 



Bae et al. 



(120091) are al so computing Bk with a mixed -action approach using HYP-smeared 
staggered valence quarks (|Hasenfratz and Knechtlil 1200 lb on the MILC ensembles. They have 
preliminary data on the coarse, fine, and superfine MILC ensembles and are computing Zg K non- 
perturbatively in the RI/MOM scheme using the Rome-Southamp ton method. When completed, 



their result should be competitive with those of RBC/UKQCD and 



Aubin et al 



(2009b) 



E. B°-B° mixing 

he mass differenc es between the heavy and light B Q , q = d,s, are given in the standard model 



by (|Buras et al. 



1990) 



f-i2 tui2 

AMf or = -^\V; q V th \\%(x t )M B jlB Bti , (206) 

where T|f is a perturbative QCD correction factor and So is the Inami-Lim function of x t = mj/M^. 
Bg q is the renormalization group invariant B® q bag parameter that can be computed in lattice QCD. 

The four-fermi operators whose matrix elements between B Q q and B q are needed to study B q 
mixing in the standard model are 

013 = [b a q a ] V -A[b c q c ]v-A , OS« = [b a q a ]s-p[b c q c }s- P , 

03l = [b a q c } S - P [b c q a } S -p, (207) 

where a,c are color indices. The leading-order B q -B q mixing matrix element is parameterized by 
the product /J : 

( g0 |OL9|fi0> MS (ju) = *-Mlf B B™{v) , (208) 
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where is related to Bg q in Eq. (|206l) in an analogous manner to Eq. (12031) . Beyond tree level, 
the operator OL q mixes with OS q , both on the lattice and in the continuum. Including the one-loop 
correction, the renormalized matrix element is given by 



(OL«) MS (v) = [l+a s -p L ^,m b )](OLT t (a)+a s -p LS (fi,m b )(OST t (a) . (209) 



,lat 



lat/ 



2M h 



20071) 



The operator 03 q is only needed to compute the width difference AT q (ILenz and Nierstd . 

The HPQCD collaboration calculated Bb C] , with q = d, s on four MILC ensembles with a ps 0. 12 
fm and two ensembl es with a . 09 fm, using an asqtad lig ht valence quark and lattice NRQCD 



for the bottom quark (IDalgic et al. , 



2007 



Gamiz et al. 



20091) . With NRQCD for the heavy quark, a 



dimension seven operator contributes to the relevant matrix el ement at order O (A® CD /Mr), which 



was also taken into account. The HPQCD collaboration finds (|Gamiz et al. 



2009) 



(210) 



fB s \JB Bs =0.266(6)(17)GeV , f Bd \jB Bi = 0.216(9) (12) GeV , 
and for the ratio 

^ = fB s ^B B ~J{fB d ^B B ~ d ) = 1.258(25)(21) , (211) 

where the errors are from statistics plus chiral extrapolation and from all other systematic errors 
added in quadrature, respectively. The chiral and continuum extrapolation is shown in the left 
panel of Fig. l37l Using the result i n Eq. (|21 1|) and the experimentally measured mass differences 



AM X , x = s,d, (lAmsler et al. 



20081) they find 

\Vtd\ 



w, 



0.214(1)(5) 



(212) 



ts\ 



where the errors are experimental and theoretical, respectively. 

A similar calcu l ation is being performed by the Fermilab Lattice and MILC collaborations 



(Evans et al. . 



2009 



20071) . They use Fermilab fermions for the heavy quarks, and, like HPQCD, 



asqtad fermions for the light valence quarks. The preliminary chiral and continuum extrapolation 
is shown in the right panel of Fig. [37] As a prelimi nary result they fin d \ = 1.205(52), with the 



statistical and systematic errors added in quadrature (jEvans et all 



2009) 



F. Hadronic contribution to the muon anomalous magnetic moment 

One of the most precisely measured quantities, and hence an astonishingly accurate test of 
QED, is the anomalous magnetic moment of the muon, a^ = (g — 2) /2. The QED contribution 
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FIG. 37 The ratio = £ ^/M B jM Bd = f Bs y / M B B B J(f Bd ^/M Bd B Bd ) as a function of the light valence quark 



collaboration (Gamiz et al. 



2009). 



mass together with rSXPT fits and the chiral and continuum extrapolation. The left panel is fro m the HPQCD 



2009) and the right panel from the Fermilab/MILC collaboration (lEvans et al. 




FIG. 38 The lowest-order diagram for the QCD correction to the m uon anomalous magnetic moment at 



0(a 2 ). The bubble represents all possible hadronic states. Figure from 



Aubin and Blum 



(120071) . 



is kno wn to four loops, with the five-loop term having been estimated — see IJegerlehnerl (|2007l 
2008|) for recent reviews. With the experimental precision to which a M is known, QCD corrections 
are important at leading order via the QCD contribution to the vacuum polarization, shown in 
Fig. m 

This leading contribution can be estimated from the experimental values o f the e + e — > hadrons 



total cross section, a^j L0 = (692.1 ±5.6) x 10 10 (IJegerlehner . 
difference between experimental and theoretical value is 



2007 



2008). Using this value the 



exp ^the 



(287 ±91) x lQ- 



ii 



(213) 
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FIG. 39 Two different rSXPT fits to U(q 2 ) for thr ee light masses: am, = 0.0031 (diamonds), 0.0062 
(squares) and 0.0124 (circles) with am s = 0.031, from 



Aubin and Blum 



(2007) which contains the details. 



about a 3.1a effect and a possible hint at effects from physics beyond the standard model. The 
leading hadronic contribution can also be estimated from x — > \ z + hadrons, giving a result of 
10 — 20 x 10~ 10 higher than from the e + e~ cross section, but this estimate is on somewhat weaker 
footing due to i so spin-breaking effects. A purely theoretical calculation of a^ L0 is thus desirable. 

The muon anomalous magnetic moment can be extracted from the full muon-photon vertex. 
The first effects from QCD, at order 0(a 2 ), ar e shown in 



vacuum polarization of the photons Tl(q 2 ) via (IBlum , 

„hlo 



2003) 



ig. [381 and can be computed from the 



Blum 



^) 2 rdtftf)ntf) 

Tt ' JO 



(214) 



(12003b . The kernel f(q 2 ) diverges as q 2 — > 0. This makes a 



with the kernel f(q 2 ) given in 
precise calculations of Tl(q 2 ) at low momentum necessary, and, in particular, makes perturbative 
co mputations unreliable. 



Aubin and Bluml (|2007l) describe such a calculation based on three MILC ensembles with lattice 
spacing a ~ 0.09 fm, and three different light quark masses. The vacuum polarization Tl(q 2 ) is 
computed from the correlator of the electromagnetic current in terms of quark fields. Aubin and 
Blum use rSXPT to fit Tl(q 2 ) at low q, (see Fig. [39b, and use the result in the integral in Eq. (12141) . 



Finally, they extrapolate to the physical light quark mass, obtaining 



,HLO 



(721 ± 15) x 10" 10 and a 



,HLO 



(748 ±21) x 10" 



10 



(215) 



132 



with a linear and quadratic fit, respectively. The errors are statistical only. Systematic errors in 
Eq. (12151) other than due to the quark mass extrapolation come from finite lattice spacing and 
finite volume effects. Given this, the lattice result should be taken as in broad agreement with the 
estimate from the e + e cross section. Further improvements need to be made before the lattice 
calculation becomes competitive with other determinations. 



G. Quark and gluon propagators in Landau gauge 



Quark and gluon propagators contain perturbative and nonperturbative information about QCD. 
Quark propagators play a crucial role in hadron spectroscopy and the study of three and four-point 
functions used in form factor and matrix element calculations. The propagators are not gauge in- 
variant, and thus have to be studied in a fixed gauge, usually the Landau gauge. Nevertheless, they 
contain gauge independent information on confinement, dynamical mass generation and sponta- 
neous chiral symmetry breaking. Quark and gluon propagators can, obviously, be studied on the 
lattice. They a re of ten treated sem i-analytically in the context of Dyson-Schwinger equations, see 



Roberts 



(2008J) and 



Fischer! (|2006r) for recent reviews. 



The Landau gauge gluo n propagator has been studied in full QCD using MILC lattices by 



Bowman et al. 



(2004. 



2007). In the continuum, the Landau gauge gluon propagator has the tensor 



structure 



(216) 



where, at tree level D(q 2 ) = l/q 2 . The bare propagator is related to the renormalized propagator 
D R (q 2 ;/u) by the renormalization condition 



D(q 2 ,a) = Z 3 (a;v)D R (q 2 ;v) , D R (q 2 ;/j)\ q2=p 2 = \ . 



(217) 



The gluon propagator in full QCD is somewhat less enhanced for momenta aro und 1 GeV than 
the quenched propagator, see Fig. |40] (left), and shows good scaling behavior (|Bowman et all 



2007|) . The gluon spectral function shows c lear vio l ations of positivity in qualitative agreement 



with Dyson-Schwinger equation studies (see 



Fischer ( 



(|2006h and references therein). 
The quark propagator has been studie d in full QCD using MILC lattice ensembles with 



lattice spacings a « 0.12 and 0.09 fm in 



Bowman et al, 



(2005b), 



Parappilly et all (|2006l) and 
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FIG. 40 The gluon dressing fu nction q 2 D(q 2 ) for quenched and dynamical configurations with a « 0.09 



fm, from 



Bowman et al. 



a « 0.12 and 0.09 fm, from 



(120071) (left), and the quar k mass function for light sea quark mass in full QCD at 



Parappillv et al\ (12006) (right). 



Furui and Nakajimal (|2006l) . The bare propagator can be parametrized, and related to the renor- 



malized propagator, by 

S(p 2 ;a) =Z(p 2 ;a)[iy-p + M(p 2 )}- 1 =Z 2 (a;iu)S R (p 2 ;iu) 



(218) 



where Z2(a;/j) = Z(p 2 ;a)\ p 2 =/J 2, and the mass function M(p 2 ) is renormalization point indepen- 



dent. Its asymptotic behavior as p — > °° is re 



condensate, see, e.g., 



Bowman et al. 



(12005al) 



ated via the OPE to the RGI quark mass and the chiral 



The quark mass function for light sea quark mass in full QCD simulations at two different lattice 
spacings is shown in Fig. |40] (right). It shows good scaling and clear indication of dynamical mass 
generation ("constituent mass") at low momenta. 



H. Further uses of MILC lattices 



Besides the calculations described in the preceding subsections, the MILC lattice ensembles 
have been used in other QCD calcul ations. These include the study of hadronic scattering lengths 



and n-body interactions, reviewed in 



Beane et al. 



(|2008al) . Furthermore, computations of nucleon 



structure, moments of parton and generalized parton distribution functions, axial nucleon cou- 
plings, electromagnetic form factors, and nuc l eon transition ampl itudes have be en done using 



MILC lattice ensembles - see 



Orginosl (12006I) . iHagled (120071) and 



Zanottil (|2009l) for recent re- 



views of lattice computations of these quantities. 
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X. FURTHER IMPROVEMENTS: A LOOK TO THE FUTURE 



While the lattice QCD simulations described in this review are quite mature, the errors of many 
of the observables computed can be reduced in various ways. Many of the calculations have omit- 
ted some of the available MILC ensembles, in particular the more challenging ones with small 
lattice spacings. Sometimes, not all the available configurations in an ensemble have been an- 
alyzed. Electromagnetic effects, where needed, have been taken from nonlattice estimates (see 
Sec.lVTl). They can be included directly in lattice simulations. Discretization effects coming from 
the fermion actions used can be further reduced by using improvements to the Fermilab action 
for heavy quarks, and by using highly improved staggered quarks for both valence and sea light 
quarks. These improvements are briefly outlined in this section. 



A. Impact of new ensembles 

The superfine (a ~ 0.06 fm) and ultrafine (a ~ 0.045 fm) ensembles listed in Table I were com- 
pleted only during the past year, as was the coarse (a ~ 0.12 fm) ensemble with three degenerate 
light quarks. The fine ensembles with mi/m s = 0.05 and with three degenerate light quarks are still 
running, but should be completed in the near future. In this paper, we have presented some prelim- 
inary results from the superfine ensembles for the hadron spectrum, the light pseudoscalar mesons 



and the topological susceptibility, and the HPQCD/UKQCD collaboration 



las recently used some 



of the superfine ensembles in its studies of charmed physics (|Daviesl . l2008h : however, the physics 
analysis of the new ensembles is in a very early stage. When it is completed, we expect these 
ensembles to have a major impact on many of the calculations described above. 

As indicated earlier, the leading finite lattice spacing artifacts for the asqtad action are of order 
a 2 /log(a). So these artifacts for the superfine and ultrafine ensembles are down from those of 
the fine ensembles by factors of 2.6 and 5.2 respectively. As one can see from Figs. Q31 [19] and 
l24l results obtained to date from the superfine ensembles are very close to the rSXPTcontinuum 
extrapolations, which should significantly reduce discretization errors in calculations that make 
use of them. Furthermore, as is illustrated in Fig.[6l the decrease in taste splitting among the pions 
with decreasing lattice spacing is consistent with a 2 / log(a) 2 , as expected. Thus, this major source 
of systematic error will be significantly reduced by use of the superfine ensembles. 

The a ~ 0.045 fm, m/ = 0.2 m s ensemble will provide an anchor point for extrapolations to 
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the continuum limit, and is particularly important for calculations which use the Fermilab method 
for heavy valence quarks. For many of these quantities the discretization errors in the heavy- 
quark action are the largest single source of systematic error. Although the size of heavy-quark 
discretization errors can be estimated using power-counting arguments, the precise form of the 
lattice spacing dependence is not explicitly known. It is thus important to have a range of lattice 
spacings in order to study the heavy-quark discretization effects. The heavy-quark errors decrease 
as a/log(a) at the worst, so we expect the 0.045 fm ensemble to reduce the heavy-quark errors 
by a factor of two in quantities of interest involving B and D mesons, which thus far have only 
been computed on ensembles with lattice spacings a ~ 0.09 fm and larger. The reduction of the 
heavy-quark discretization errors does not require the full set of light-quark masses that we have 
calculated at coarser lattice spacings; thus, we have generated only one ensemble at a ~ 0.045 fm. 
By including the superfine and ultrafine ensembles into our work on heavy-light mesons, in con- 
junction with improving the statistics, we expect to determine the leptonic decay constants, the 
mixing parameters and the corresponding semileptonic form factors to an accuracy of better than 
5%. 

The physical strange quark mass is not light enough for chiral perturbation theory to converge 
rapidly in its vicinity. To anchor chiral fits and to test the convergence of chiral perturbation theory, 
it is therefore extremely helpful to have ensembles with the strange sea quark mass held fixed at 
a value well below the physical strange quark mass. Furthermore, with three dynamical quark 
flavors, there are two interesting chiral limits to be considered: the two-flavor limit, in which the 
u and d quarks become massless while the s stays at its physical mass, and the three-flavor chiral 
limit, where all three quarks become massless. The difference of various quantities in these two 
limits is an important probe of the nature of chiral symmetry breaking in QCD. The extrapolation 
to m s = necessary for the three-flavor chiral limit is a long one, with attendant large errors. 
The new ensembles with three degenerate light quarks were created to help address these issues. 
We estimate that incorporating all the superfine ensembles into the analysis, as well as all the 
configurations with the strange sea quark mass held fixed below its physical value, will allow us 
to reduce the systematic errors on f % and fx to 2% or better, and should dramatically reduce the 
errors in low energy constants and quantities such as the ratio of the two flavor to three flavor 
condensates, (uu)2/ (mm) 3. This would be an important milestone for lattice QCD calculations. We 
also expect corresponding improvements in other physical quantities of interest. In particular, our 
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evaluation of \V US \ should become significantly more accurate than the current world average. 



B. Electromagnetic and isospin breaking effects 



Most lattice calculations have not included electromagnetic or isospin breaking effects. How- 
ever, as the precision of calculations increases, including these effects will become increasingly 
important. In fact, we have already seen in Sec. [VI] that electromagnetic effects are important in 
the determination of the u and d quark masses. Another interesting challenge for lattice QCD 
would be to determine the proton-neutron mass difference, which will require accounting for the 
differences of both the u and d quark masses and their charg es. 



The pioneering work by 



Duncan. Eichten. and Thacken ( 



1996 . 



1997j) regarding electromag- 



netic effects was done with quenched U(l) and quenched SU(3) fields. More recently, the RBC 
collaboration has be en p ursuing such calcu lations but with domain-wall dynamical quarks. In 

(12007b. electromagnetic effects on 7t and K meson masses 



Yamada et al. 



(2006J) and 



Blum et al 



were calculated in Nf = 2 configurations. 



Beane. Orginos. and Savage! (|2007bh have used MILC 



configurations with a « 0.12 fm to study isospin breaking for the nucleons using domain-wall 
valence quarks. 

Electroma gnetic effects i n lowest order chiral perturbation theory were first studied some 40 
years ago by iDashenl (| 19691) . A key result known as Dashen's Theorem is that electromagnetic 
splittings of the pions and kaons are equal at this order, i.e., 



AM 2 = AM 2 - AM 2 = (M| ± - M| ) em - (mJ ± - M 2 ) , 



(219) 



vanishes. 



Recently, iBijnens and Danielssonl (|2007h have calculated electromagnetic corrections in par- 
tially quenched perturbation theory, which are particularly pertinent for analysis of lattice QCD 
calculations. They have emphasized that a combination of meson masses with varying charges and 
quark masses is a very close approximation to AMj^: 

AM 2 = M 2 (xi,X3,4i,43)-M 2 (xi,X3,43,4 3 ) 

- ^ 2 (%i,%i^i^3)+M 2 (Xi,%i^3,^)- (220) 



Here %, = 2Bm qi , where B is the continuum version of the low energy constant defined in Eq. (1391) , 
and qi is the quark charge. In their notation, i = 1 (3) refers to the valence u (d) quark, respectively. 
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FIG. 41 CoiTection to Dashen's theore m, as a func t ion of the LO 71 mass squared (equivalent to the pion 
mass squared with e 2 = 0). Figure from] 



Basak et al. 



(2009) 



MILC ha s rece ntly begun to explore electromagnetic effects on the pseudoscalar masses 



(|Basak et al. . 



2009), using the quenched approximation for electromagnetism. The initial study 
on a ~ 0.15 fm ensembles yielded promising results. The key result is a rough estimate of the 
correction to Dashen's theorem. In Fig. EH we show results for two dynamical ensembles for 
various light valence masses. After fitting the results and performing the chiral extrapolation, 
we find that 0.7 x 10~ 3 GeV 2 < AMp < 1.8 x 10 ~ 3 GeV 2 . A recent phenomenological estimate is 



1.07 x 10 3 GeV (|Bijnens and Danielssonl 



2007). 



It will be very interesting to extend this work to smaller lattice spacings and eventually to 
include dynamical electromagnetic effects. There is also the prospect of including isospin breaking 
in the generation of the configurations. 



C. Heavy Wilson fermion improvement program 



The leading discre tization errors contained in the Wilson/clover action applied to heavy quarks 
have been analyzed in lOktay and Kronfeldl (12008b . in an extension to the original Fermilab formal- 
ism. Since the heavy quarks introduce an additional scale 1/mg, they consider all the operators 
which have power counting of X 3 (k ~ Aa or A/mq) and v 6 for the heavy-light (HQET) and 
heavy-heavy (NRQCD) systems, respectively. This leads to actions containing all possible dimen- 
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sion six and some dimension seven operators. Many of these are redundant and may be chosen 
for calculational convenience by considering field transformations. For example, multihop time 
derivative operators (which spoil nice properties of the transfer matrix) may be eliminated in this 
way. Tree-level matching of observables in the continuum and lattice QCD actions shows that six 
new operators beyond the original Fermilab action are required at this level of improvement, four 
of dimension six and two of dimension seven. In all, there are a total of nineteen nonredundant 
operators at this level, and one-loop matching will presumably introduce more of these. One can 
estimate the uncertainties due to nonzero lattice spacing by calculating the mismatch between the 
lattice short-distance coefficients and their continuum counterparts. Initial estimates show that the 
new lattice action reduces the errors to the few-percent level. 



D. Preliminary studies of the HISQ action 



As discussed in Sec. [HI the HISQ action improves taste symmetry and is well suited for future 
studies with dynamical quarks. Subtleties with dynamical HISQ simulations, in particular from 
th e reunitarizatiqn step ? Eq. (f85l) . which can lead to large contributions to the force, are described 



in 



Bazavov et al. 



(l2009h 



The first study of how the HISQ action reduces the splitting between different ta stes of pions 



was undertaken by the HPQCD and UKQCD collaborations in 



Follana et al. 



(120071) . They used 



valence HISQ on the asqtad sea quark configura t ions g enerated by MILC. Similar findings for 



HISQ sea quarks were reported in 



Bazavov et al. 



(|2009|) . The results of a more recent study are 



summarized in Fig .1421 The splittings between the Goldstone and the other pion tastes for the HISQ 
action are reduced by a factor of 2.5-3 compared to asqtad (notice a vertical line that indicates a 
factor of 3 in logarithmic scale in Fig. [42b. Two HISQ ensemble s, with a 0.09 and .12 fm, are 



shown. The difference between the results presented here and in 



Bazavov et al. 



( 20091) is that the 



current study uses t he improved gauge action with the one-loop fermion corrections induced by 



the HISQ fermions (|Hart et all 



2009a. b), and the ensembles were tuned to be close to the line of 



constant physics with mi = 0.2m s . 
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FIG. 42 The taste splittings as function of a 2 a 2 for the asqtad and HISQ actions (with the latter indicated 
by dashed boxes). 

XL SUMMARY AND CONCLUSIONS 

There has been a dramatic improvement in the accuracy of lattice QCD calculations over the 
past decade due to a combination of developments: 

• The use of improved actions significantly reduces finite lattice spacing artifacts, greatly im- 
proving the accuracy of extrapolations to the continuum limit. The asqtad improved stag- 
gered quark action the MILC collaboration has used provides a particularly strong reduction 
in taste symmetry breaking, the most challenging finite lattice spacing artifact for staggered 
quarks. The HISQ action improves on asqtad in this respect by an additional factor of three. 
In general, one finds that a HISQ ensemble has lattice artifacts approximately half the size 
of an asqtad ensemble with the same lattice spacing. 

• The inclusion of up, down and strange sea quarks with realistic masses is critical for reducing 
errors to the few percent level, as is illustrated in Fig. [TJ 

• The use of partially quenched chiral perturbation theory and, for staggered quarks, rooted 
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staggered chiral perturbation theory have greatly improved the accuracy of the extrapolation 
of lattice data to the physical masses of the up, down and strange quarks. 

• Improved algorithms, such as RHMC, have enabled the generation of gauge field ensembles 
with significantly smaller lattice spacings and lighter quark masses than had previously been 
possible. These new algorithms have changed the balance between gauge field configuration 
generation and physics analysis on the configurations. Whereas the former used to take the 
bulk of the computing resources, now the resources required for an analysis project often 
rival those that went into the generation of the configurations. 

• The vastly increased computing resources available to lattice gauge theorists over the past 
decade have enabled us to take advantage of the developments enumerated above. For exam- 
ple, between 1999 and 2008, the total floating point operations used per year by the MILC 
Collaboration increased by approximately three orders of magnitude. 

The MILC collaboration has taken advantage of these developments to generate, over the past 
ten years, the ensembles of asqtad gauge field configurations detailed in Table IB This is the first 
set of ensembles to have a wide enough range of small lattice spacings and light quark masses to 
enable controlled extrapolations of physical quantities to the continuum and chiral limits. These 
ensembles are publicly available, and we and others are using them to calculate a wide range of 
physical quantities of interest in high energy and nuclear physics. This work has included calcula- 
tions of the strong coupling constant, the masses of light quarks and hadrons, the properties of light 
pseudoscalar mesons, the topological susceptibility, the masses, decays and mixings of heavy-light 
mesons, the charmonium and bottomonium spectra, the ^° — K° mixing parameter Bk, the mass 
of the B c meson, the it — n and N — N scattering lengths, generalized parton distributions, and 
hadronic contributions to the muon anomalous magnetic moment. The errors in these quantities 
have typically decreased by an order of magnitude as the library of ensembles has grown, with 
further improvements expected as the superfine and ultrafine ensembles are fully analyzed, and 
HISQ ensembles become available. 

A number of quantities have been calculated to an accuracy of a few percent, and some predic- 
tions have been made that were later verified by experiment. The work of the Fermilab Lattice, 
MILC and HPQCD/UKQCD collaborations on the decays and mixings of heavy-light mesons and 
the decays of light pseudoscalar mesons has reached a level of accuracy where it is having a signifi- 
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cant impact on tests of the standard model and the search for new physics. However, high precision 
has been obtained only for quantities that are most straightforward to calculate. There are many 
quantities, such as scattering phase shifts, the masses and widths of hadrons that are unstable under 
the strong interactions, and parton distribution functions, which are of great interest, but continue 
to pose major challenges. 

Because it is relatively inexpensive to simulate, the asqtad quark action was the first to produce 
a set of gauge field ensembles with a wide enough range of lattice spacings and sea quark masses 
to enable controlled extrapolations to the continuum and chiral limit. However, such ensembles are 
also being produced with other quark actions, such as Wilson-clover, twisted mass, domain wall 
and overlap. These ensembles are already producing impressive results. Over the next few years 
one can expect major advances on a wide variety of calculations with critical checks coming from 
the use of different lattice formulations of QCD. Finally, the techniques that have been developed 
for the study of QCD can be applied to study many of the theories that have been proposed for 
physics beyond the standard model. Such work is just beginning, but appears to have a very bright 
future. 
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